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ABSTRACT 


We  develop  models  and  solution  methodologies  to  solve  the  discrete-time  path- 
optimization  problem  where  a  single  or  multiple  searchers  look  for  a  moving  target 
in  a  hnite  set  of  cells.  The  single  searcher  is  constrained  by  maximum  limits  on 
the  consumption  of  several  resources  such  as  time,  fuel,  and  risk  along  any  path. 
We  develop  a  specialized  branch-an-bound  algorithm  for  this  problem  that  utilizes 
several  new  network  reduction  procedures  as  well  as  new  bounding  technique  based 
on  Lagrangian  relaxation  and  network  expansion.  The  resulting  algorithm  is  quite 
efficient  and  promising.  For  the  multiple  searchers,  an  optimal  set  of  paths  (search 
plan)  is  determined  by  taking  advantage  of  the  cooperative  search  effect.  We  present 
a  new  exact  algorithm  and  two  new  heuristics  to  hnd  an  optimal  or  near-optimal 
search  plan.  One  of  the  heuristics  is  based  on  the  cross-entropy  method  and  is  found 
to  perform  well  for  a  broad  range  of  problem  instances.  The  exact  algorithm  deals  with 
the  specihc  case  of  homogeneous  searchers  and  is  based  on  outer  approximations  by 
several  new  cutting  planes.  In  addition,  we  prove  that  under  certain  assumptions  the 
path-optimization  problem  becomes  equivalent  to  a  large-scale  linear  mixed-integer 
program. 
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EXECUTIVE  SUMMARY 


This  dissertation  develops  models  and  solution  methodologies  to  solve  the  dis¬ 
crete  time-and-space  path-optimization  problems  where  a  single  searcher  or  multiple 
searchers  look  for  a  non-evading  moving  target.  We  extended  the  single  searcher 
problem  (SSP),  generally  called  the  optimal  searcher  path  problem,  to  realistic  situ¬ 
ation  where:  (1)  the  searcher  needs  to  change  its  flight  altitude  during  a  mission  to 
balance  sensor  quality  and  risk  from  ground  threats;  (2)  the  searcher  is  subject  to 
constraints  on  multiple  resources  such  as  fuel  consumption  and  risk  exposure  along 
any  path;  and  (3)  the  search  effect  depends  on  the  current  and  previous  locations  of 
the  searcher.  The  latter  situation  may  occur  if  moving  from  some  search  sector  to 
a  new  sector  results  in  a  lower  search  effect  due  to  ineffective  search  during  transit. 
We  refer  to  this  extended  SSP  as  the  resource-constrained  single  searcher  problem 
(RSSP). 

We  present  new  bounding  techniques  (static  bound  and  directional  static 
bound)  and  network  reduction  procedure  in  a  branch-and-bound  algorithm  to  opti¬ 
mally  solve  the  single  searcher  problem  (SSP).  The  static  bound  simplihes  the  bound 
calculation  substantially  at  the  expense  of  a  weaker  bound.  The  directional  static 
bound  improves  the  weaker  static  bound.  The  network  reduction  procedure  takes 
advantage  of  the  initial  positions  of  searcher  and  target,  and  eliminates  parts  of  the 
branch-and-bound  tree  while  retaining  a  optimal  solution.  The  resulting  algorithm 
solves  a  problem  instance  with  15  by  15  cells  and  a  time  horizon  of  20  optimally  in 
less  than  4  seconds  on  average,  which  is  a  signihcant  reduction  from  the  9  minutes 
required  by  a  state-of-the-art  algorithm. 

We  extend  our  branch-and-bound  algorithm  for  the  SSP  to  treat  the  RSSP 
and  construct  a  new  bounding  technique  based  on  Lagrangian  relaxation.  In  addition, 
we  develop  novel  network  reduction  techniques  using  dominance  tests.  The  resulting 
algorithm  solves  a  realistic  problem  instance  (10  by  10  cells,  two  altitudes,  time 


xvn 


horizon  40,  and  two  resource  constraints),  on  average,  in  less  than  10  minutes. 

We  extend  the  specialized  branch-and-bound  algorithm  to  the  case  of  multi¬ 
ple  searchers  and  also  develop  two  new  heuristics  (static  bound  heuristic  and  cross¬ 
entropy  heuristic).  In  the  context  of  search  problems,  this  dissertation  appears  to  be 
the  first  one  which  utilizes  the  cross-entropy  method.  Among  these  three  algorithms, 
the  cross-entropy  heuristic  performs  best  for  a  broad  range  of  problem  instances. 
While  the  branch-and-bound  algorithm  and  the  static  bound  heuristic  are  limited 
to  small  problem  instances,  the  cross-entropy  heuristic  obtains  good  approximate 
solutions  of  moderately  sized  instances  (15  by  15  cells,  time  horizon  18,  and  three 
searchers)  in  about  10  minutes. 

For  the  case  with  multiple  identical  searchers,  we  construct  an  exact  algorithm 
based  on  outer  approximations  by  cutting  planes.  We  introduce  several  new  cutting 
planes  (multiple-cut,  strong-cut,  relative-cut,  and  symmetric-cut),  combine  them,  and 
develop  a  specialized  outer  approximation  algorithm.  For  a  moderate-size  problem 
instance  (15  by  15  cells,  time  horizon  16,  and  three  searchers),  the  resulting  algorithm 
essentially  provides  an  approximate  solution  that  is  guaranteed  to  be  within  5%  of  an 
optimal  solution  in  less  than  about  20  minutes.  Instances  with  more  searchers  (such 
as  10,  15  and  30  searchers)  are  solved  even  faster.  In  addition,  we  prove  that  under 
certain  assumptions  the  nonlinear  convex  multiple  homogeneous  searcher  problem  is 
equivalent  to  a  large-scale  linear  mixed-integer  program. 


I.  INTRODUCTION 

A.  MOTIVATION 

The  need  to  search  for  moving  objects  arises  in  many  civilian  and  military 
applications.  Rescue  teams  search  for  lost  persons  and  shipwrecks.  Autonomous 
robots  could  hnd  victims  and  survivors  after  natural  disasters.  In  military  search 
operations,  patrol  aircraft  and  unmanned  aerial  vehicles  (UAVs)  look  for  high-valued 
targets  such  as  missile  launchers  and  terrorist  vehicles.  In  all  these  applications, 
the  choice  of  paths  for  the  searchers  strongly  influences  the  probability  of  hnding 
the  targets  within  a  specihc  time  limit.  Unfortunately,  the  problem  of  selecting  the 
“best”  path  is  fundamentally  hard  due  to  the  nonlinearity  induced  by  the  probability 
of  detection.  For  example,  looking  twice  by  a  searcher  generally  does  not  double  the 
detection  probability. 

In  military  search  missions  with  UAVs,  the  problem  of  hnding  good  paths 
for  searchers  (UAVs)  is  complicated  even  further.  Rather  than  manned  aircraft, 
UAVs  are  preferred  to  operate  in  dangerous  environments  such  as  hostile  regions  and 
battlehelds.  In  these  applications,  the  searchers  (UAVs)  often  need  to  balance  search 
effectiveness  with  threats  to  themselves  posed  by  surface-to-air  missiles  and  small- 
arms  hre.  Moreover,  the  search  problem  becomes  signihcantly  more  difficult  as  the 
number  of  searchers  grows.  With  the  technological  advances  in  cooperative  control  of 
multiple  UAVs,  such  multi-searcher  problems  are  appearing  with  increasing  frequency 
in  applications.  The  need  for  optimal  path  planning  in  these  complicated  situations 
motivates  this  dissertation. 

In  this  dissertation,  we  consider  several  search  problems  including  situations 
with  searchers  subject  to  threats  as  well  as  multiple  searchers.  We  particularly  focus 
on  search  by  aerial  assets  that  are  subject  to  threats  from  the  ground.  We  formulate 
search  problems  for  single  and  multiple  searchers  and  develop  solution  methodologies 
for  efficiently  hnding  an  optimal  or  near-optimal  path. 
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B.  SCOPE  OF  DISSERTATION 

We  consider  discrete  time-and-space  path-optimization  problems  where  a  non¬ 
evading  target  moves  in  a  two-dimensional  area  of  interest  (AOI).  The  AOI  is  dis¬ 
cretized  into  a  hnite  set  of  cells  and  the  target  moves  between  the  cells  in  discrete 
time.  Stone  [47]  and  Brown  [7]  introduced  the  basic  forms  of  this  search  problem. 
A  single  searcher  or  multiple  searchers  move,  in  discrete  time,  through  a  discretized 
two-dimensional  space  (on  the  ground  in  the  AOI  or  at  a  constant  altitude  over  the 
AOI)  or  a  discretized  three-dimensional  airspace  over  the  AOI.  The  searchers’  space 
is  represented  by  a  directed  graph,  where  vertices  represent  waypoints  (or  search 
sectors)  from  which  a  searcher  can  search  a  particular  cell  in  the  AOI  and  where  di¬ 
rected  edges  represent  transition  between  two  waypoints.  Thus  a  path  of  the  searcher 
is  represented  as  a  sequence  of  waypoints  (vertices). 

We  note  that  our  aim  is  not  to  hnd  an  optimal  flyable  trajectory  that  fully 
accounts  for  turn  radius,  climb/dive  angle,  and  max/min  speed  for  the  searcher.  We 
aim  for  “high-level”  control  of  the  searchers  where  waypoints  may  represent  points 
in  space  that  are  many  minutes  apart  in  flying  time  and  may  represent  areas  to  be 
searched  for  some  period  of  time.  We  consider  imperfect  searchers  that  may  search  a 
cell  that  contains  a  target  without  hnding  the  target.  However,  in  this  dissertation, 
we  assume  that  the  searchers  will  not  report  a  target  in  a  cell  where  there  is  no  target. 

The  goal  for  the  searchers  is  to  maximize  the  probability  of  detecting  the  target 
within  a  hnite  mission  time  or,  equivalently,  to  minimize  the  probability  that  the 
target  is  not  detected  during  the  duration  of  the  mission.  The  searchers  are  subject 
to  constraints  on  path  continuity,  mission  time,  and  possibly  other  “resources”  such 
as  risk  exposure  and  fuel  consumption. 

This  dissertation  refers  to  the  path  optimization  problem  for  a  single  searcher 
with  path-  and  time-constraints  as  the  single  searcher  problem  (SSP).  Other  authors 
refer  to  this  problem  as  the  optimal  searcher  path  problem  [41,  49,  32].  In  this 
dissertation,  we  focus  on  the  following  variants  of  SSP. 
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1.  SSP  with  additional  constraints  on,  e.g.,  risk  exposure  to  threats  and  fuel 
consumption  during  the  missions.  We  refer  to  this  path-optimization  problem 
as  the  resource-constrained  single  searcher  problem  (RSSP). 

2.  SSP  for  multiple  searchers  where  the  probability  of  detecting  the  target  by  at 
least  one  searcher  is  maximized.  We  term  this  problem  the  multiple-searcher 
problem  (MSP). 

3.  MSP  for  homogeneous  searchers  where  the  fact  that  searchers  are  identical  is 
utilized.  We  refer  to  this  multiple-homogeneous-searcher  problem  as  MHSP. 

This  dissertation  always  assumes  that  the  searchers  are  subject  to  constraints  on 
path  continuity  and  mission  time.  We  note  that  the  multiple-searcher  problem  with 
resource-constraints  (risk  exposure  and  fuel  consumption,  etc.)  is  outside  of  the  scope 
of  this  dissertation. 

C.  LITERATURE  SURVEY 

The  topic  of  searching  for  a  moving  target  in  discrete  time  and  space  where 
a  searcher’s  path  is  constrained  has  received  much  attention.  Benkoski  et  al.  [41] 
provide  a  comprehensive  survey  of  this  problem  (until  1990).  The  path-  and  time- 
constrained  search  problem  with  a  stationary  target  is  NP-complete  [49].  Thus  the 
path-  and  time-constrained,  moving-target,  search  problem  is  at  least  as  difficult. 

One  scheme  to  hud  an  optimal  path  for  a  searcher  is  the  branch-and-bound 
procedure.  A  branch-and-bound  procedure  was  hrst  studied  by  Stewart  [46].  Since 
Stewart’s  “bounds”  are  not  valid,  an  optimal  branch  of  the  enumeration  tree  may  be 
fathomed.  More  recent  studies  [18,  34,  52,  32]  focus  on  the  development  of  specialized 
branch-and-bound  algorithms  for  hnding  an  optimal  path  for  the  single  searcher.  In 
these  algorithms,  a  path  is  a  sequence  of  vertices  that  the  searcher  will  visit  and 
branching  corresponds  to  extending  a  subpath  by  one  more  vertex.  Bounds  on  the 
optimal  value  of  the  problem  are  obtained  by  replacing  the  probability  of  detection 
with,  effectively,  the  expected  number  of  detections  [34,  52,  32].  Bounds  are  also 
obtained  by  assuming  that  the  searcher  can  divide  its  effort  among  multiple  cells  each 
time  period  [18].  The  efficiency  of  various  branch-and-bound  methods  is  investigated 
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in  [52] ,  with  emphasis  on  the  tradeoff  between  the  accuracy  of  the  bound  employed 
and  the  time  required  to  compute  it. 

An  alternative  approach  to  obtain  an  optimal  path  is  based  on  dynamic  pro¬ 
gramming  and  partially  observable  Markov  decision  processes  [17].  The  approach 
can  solve  small  problems  quickly,  but  requires  a  large  amount  of  computer  storage  as 
the  problem  size  increases.  Later  the  same  author  reports  that  a  branch-and-bound 
algorithm  [18]  performs  more  efficiently  and  can  solve  larger  problems  than  the  exist¬ 
ing  dynamic  programming  procedures.  Thus,  we  mainly  focus  on  branch-and-bound 
procedures  in  this  dissertation. 

1.  Resource-constrained  Single  Searcher  Problem 

In  the  studies  listed  above,  it  is  assumed  that  the  searcher  moves  on  the 
ground  in  the  AOI  or  searches  the  AOI  from  a  constant  altitude.  However,  in  mil¬ 
itary  search,  surveillance,  and  reconnaissance  operations  over  land,  factors  such  as 
risk  and  fuel  consumption  become  independent  elements  of  concern  for  planners  [37], 
and  the  searcher  is  required  to  change  its  flight  altitude  during  a  mission.  Risk  to 
the  searcher  arises  from  exposure  to  enemy  threats  on  the  ground  such  as  small-arms 
hre,  anti-aircraft  artillery,  and  surface-to-air  missiles.  Risk  of  pilot  error  (e.g.,  for  a 
low-flying  helicopter  [29])  and  mechanical  failure  (e.g.,  for  small  and  unreliable  UAVs 
[31])  may  also  be  significant.  Fuel  consumption  tends  to  be  proportional  to  the  search 
duration.  However,  for  small  UAVs,  it  is  not  always  proportional  to  the  mission  time, 
due  to  variation  in  the  searcher’s  speed  and/or  altitude,  as  well  as  varying  weather 
conditions.  During  search  over  land,  terrain  features  require  altitude  changes.  The 
altitude  is  also  varied  to  balance  the  searcher’s  risk  with  image  quality.  We  note 
that  image  quality  is  of  particular  concern  for  small  UAVs  operating  with  low-  to 
moderate-quality  sensors  (For  information  about  UAV  surveillance  operations  and 
sensors,  see  [37,  42]).  In  civilian  search  and  surveillance  operations  over  land,  we  may 
face  many  of  the  same  factors  with  the  exception  of  ground  fire. 


4 


There  are  only  a  few  studies  dealing  with  constraints  on  risk  and  fuel  con¬ 
sumption  for  the  searchers.  Thomas  and  Eagle  [48]  consider  a  search  problem  with 
a  random  time  horizon;  i.e.,  with  a  perishable  target  with  a  random  lifetime  that 
follows  a  particular  distribution  (a  geometric  distribution  in  that  study).  This  situa¬ 
tion  is  similar  to  having  the  searcher  being  subject  to  threats  and  having  the  search 
terminate  after  searcher  “failure.”  We  will  consider  a  similar  situation,  but  will  limit 
the  total  risk  a  searcher  can  be  exposed  to  during  the  mission.  Secrest  [43]  considers 
optimal  path  planning  with  multiple-objectives  (simultaneously  accounting  for  flight 
time,  risk  and  fuel  consumption,  etc.)  in  surveillance  missions  while  visiting  all  as¬ 
signed  locations.  The  resulting  model  does  not  consider  the  probability  of  detecting 
a  target. 

Models  for  military  aircraft  routing  frequently  consider  minimum  risk  along  a 
path  subject  to  fuel  consumption  and  mission  duration  constraints,  see,  e.g.,  [35,  55, 
10].  These  models  have  similar  constraints  to  the  ones  in  RSSP,  but  deal  with  linear 
and  time-invariant  objective  functions.  In  contrast,  RSSP  has  a  nonlinear,  time- 
variant  objective  function  representing  the  probability  of  detection  (see  Chapter  II). 
Only  recently  have  researchers  extended  the  models  for  aircraft  routing  to  situations 
with  nonlinear  objective  functions,  and  then  only  for  special  cases  of  path-dependent 
risks  [29]. 

2.  Multiple-Searcher  Problem 

The  branch-and-bound  procedure  is  also  applicable  in  solving  the  multiple- 
searcher  problems  (MSP/MHSP).  Dell  et  ah  [13]  present  an  exact  procedure  (branch- 
and-bound  algorithm)  and  six  heuristics  (local  search,  expected  detection  heuristics, 
genetic  algorithms,  and  moving  time-horizon  heuristic)  to  solve  the  MSP/MHSP. 
However,  the  branch-and-bound  algorithms  need  prohibitive  runtime  to  guarantee  an 
optimal  path  for  multiple-searchers,  because  the  number  of  possible  paths  grows  expo¬ 
nentially  with  the  number  of  searchers.  This  motivates  the  development  of  heuristic 
algorithms.  We  find  other  heuristic  algorithms  based  on  myopic  and  receding-horizon 
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approaches  in  [23]  and  sequential  optimization  of  each  searcher  in  [26] .  For  a  specihc 
scenario,  Riehl  et  ah  [28]  present  a  receding-horizon  approach  that  jointly  optimizes 
paths  and  sensor  orientations  for  multiple  searchers.  For  a  stationary-target  case, 
Song  et  ah  [45]  use  a  forward  induction  approach  to  hnd  an  optimal  search  strategy. 

Control  and  robotic  communities  have  analyzed  the  search  problems  for  mul¬ 
tiple  UAVs  and  sensors  as  well.  Ryan  et  al.  [2]  review  the  current  issues  and  develop¬ 
ments  in  the  held  of  cooperative  control  (real-time  navigation,  collision  avoidance,  and 
hight  formation,  etc.).  Decentralized  search  techniques  [20,  54]  and  Bayesian  heuristic 
approaches  [53,  25]  are  also  studied  in  the  context  of  real-time  search  operations. 

D.  CONTRIBUTIONS 

This  study  considers  the  path  optimization  problems,  including  “searcher  sub¬ 
ject  to  threat”  and  cases  with  multiple  searchers.  We  extend  the  classical  single 
searcher  problem  (SSP)  to  realistic  situations  such  as:  (1)  the  searcher  needs  to 
change  its  hight  altitude  during  a  mission  to  balance  sensor  quality  and  risk  from  the 
ground  treats;  (2)  the  searcher  is  subject  to  constraints  on  multiple  resources  such  as 
fuel  consumption  and  risk  exposure  along  its  path;  and  (3)  the  search  ehect  depends 
on  the  current  and  previous  locations  of  the  searcher.  The  latter  situation  may  occur 
if  moving  from  some  search  sector  to  a  new  sector  results  in  a  lower  search  ehect  due 
to  inehective  search  during  transit. 

We  present  new  bounding  techniques  and  network  reduction  procedures  in  a 
branch-and-bound  algorithm  to  efficiently  solve  the  SSP.  Based  on  this  algorithm,  we 
develop  a  specialized  branch-and-bound  procedure  to  solve  the  resource-constrained 
single  searcher  problem  (RSSP)  by  utilizing  several  novel  network  reduction  proce¬ 
dures  as  well  as  new  bounding  technique  taking  account  of  the  multiple  resource 
constraints.  We  also  develop  several  new  algorithms  to  solve  the  multiple  searchers 
problems  (MSP/MHSP). 
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E.  ORGANIZATION 

The  reminder  of  this  dissertation  is  outlined  as  follows.  In  the  next  chapter, 
we  formulate  RSSP  and  develop  a  specialized  branch-and-bound  algorithm  for  RSSP. 
The  resulting  algorithm  utilizes  a  new  bounding  technique,  applies  several  network 
reduction  procedures,  and  exhibits  promising  behavior  in  computational  tests  on  in¬ 
stances  of  SSP  and  RSSP.  In  Chapter  III,  we  formulate  MSP  by  extending  SSP  and 
present  three  new  algorithms  (an  exact  procedure  and  two  heuristics)  to  solve  MSP 
and  examine  their  computational  efficiency.  In  Chapter  IV,  we  focus  on  MHSP.  We 
use  the  convexity  of  the  nonlinear  objective  function  (the  non-detection  probability) 
and  construct  an  exact  algorithm  using  cutting  planes  (outer  approximations).  We 
prove  that  under  certain  assumptions  the  problem  is  equivalent  to  a  large-scale  lin¬ 
ear  mixed-integer  program.  We  also  introduce  several  new  cuts  for  the  MHSP  and 
demonstrate  their  effect  in  computational  tests.  Chapter  V  provides  conclusions  of 
our  study  as  well  as  suggestions  for  future  work. 
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II.  RESOURCE-CONSTRAINED  SINGLE 
SEARCHER  PROBLEM 

This  chapter  formulates  the  resource-constrained  single  searcher  problem  (RSSP) 
by  generalizing  existing  models  of  the  single  searcher  problem  (SSP)  to  the  case 
with  multiple  resource  constraints.  The  RSSP  also  considers  an  altitude-dependent 
searcher  with  its  search  effect  depending  on  not  only  the  searcher’s  current  loca¬ 
tion  but  also  its  previous  location.  We  combine  the  branch-and-bound  methodology 
for  solving  SSP  with  the  line  of  research  on  the  resource-constrained,  shortest-path 
problem  [10].  Specihcally,  we  merge  the  algorithms  in  [32]  and  [10],  and  develop  a  spe¬ 
cialized  branch-and-bound  algorithm  for  RSSP.  The  resulting  algorithm  utilizes  new 
bounding  techniques,  applies  several  new  network  reduction  procedures,  and  exhibits 
promising  behavior  in  computational  tests  on  instances  of  SSP  and  RSSP. 

A.  PROBLEM  DESCRIPTION  AND  FORMULATION 

The  area  of  interest  (AOI)  is  discretized  into  a  hnite  set  of  cells  C  =  {1, . . . ,  C*} 
(see  Figure  1)  and  the  time  horizon  is  discretized  into  a  hnite  set  of  time  periods 
T  =  {1,2,  ...,T}.  A  target  occupies  one  cell  each  time  period  and  moves  among 
cells  according  to  a  Markov  process  with  known  transition  matrix  P  =  (Pc,c'))  where 
Pc,c'  is  the  probability  that  the  target  moves  from  cell  c  to  cell  c'  during  one  time 
period.  This  Markovian  target  motion  is  hrst  modeled  in  Stone  [47].  Let  p{-,t)  = 
[p{l,t),p{2,t), . . .  ,p{C,t)],  where  p{c,t)  is  the  probability  that  the  target  is  in  cell 
c  G  C  at  the  beginning  of  time  period  t  E  T  and  the  target  has  not  been  detected 
before  t.  We  refer  to  p{-,t)  as  the  undetected  target  distribution.  We  note  that  p{-,t) 
may  not  be  a  probability  mass  function  as  the  sum  of  the  elements  could  be  less  than 
1.  The  initial  target  distribution  p{-,  1)  is  known. 

A  single  searcher  moves  through  a  designated  airspace  over  the  area  of  interest 
with  the  goal  of  hnding  the  moving  target  on  the  ground.  The  airspace  over  each 
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Figure  1.  A  discretized  area  of  interest  composing  of  C*  =  4  cells. 

cell  c  G  C  is  vertically  discretized  into  a  set  of  altitudes  (or  boxes)  Td  =  {1,  2, . . . ,  H} 
(see  Figure  2).  For  any  c  G  C  and  h  eH,  we  refer  to  the  cell-altitude  pair  (c,  h)  as  a 
waypoint  where  the  searcher  can  loiter  and  carry  out  search  of  cell  c.  We  assume  that, 
from  waypoint  (c,  h) ,  the  searcher  can  sweep  cell  c  only.  We  model  the  designated 
airspace  by  a  directed  network  (V,£^),  with  set  of  vertices  V  and  set  of  directed 
edges  S,  in  which  vertices  v  =  {c,h)  E  V  represent  waypoints  and  directed  edges 
e  =  (u,  v')  E  S  represent  transition  between  waypoints  v,  v'  E  V.  The  waypoints 
v,v'  eV  are  adjacent  to  each  other  if  v  and  v'  touch  at  any  portions.  The  searcher 
can  only  transit  between  two  waypoints  that  are  adjacent  to  each  other.  Let  J^{v)  eV 
be  the  set  of  vertices  that  are  adjacent  to  u  G  V.  We  refer  to  J^{v)  as  the  forward 
star  of  vertex  v.  We  adopt  the  convention  that  v  E  J^{v)  for  all  v  eV.  Then,  the  set 
of  edges  £  =  {(u,u')  G  V  x  V|  F  G 

During  each  time  period  t  E  T,  the  searcher  is  at  a  particular  vertex  (way- 
point).  We  assume  there  is  no  transit  time  between  waypoints.  Hence,  {v,v')  E  £ 
simply  represents  search  at  waypoint  v  followed  by  search  at  waypoint  v'  in  the  next 
time  period.  The  situation  with  nonzero  transit  time  between  waypoints  can  be  mod¬ 
eled,  at  least  approximately,  by  introducing  artificial  vertices.  (We  refer  to  [32]  for 
a  comprehensive  study  of  nonzero  transit  times.)  We  note  that  the  edge  {v,v)  E  £ 
represents  searching  at  waypoint  v  for  two  consecutive  time  periods. 
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Figure  2.  A  discretized  airspace  over  the  area  of  interest  (Figure  1)  with  H  =  2 
altitudes. 

Let  0  :  V  — >■  C  be  the  function  that  specifies  the  cell  over  which  a  vertex  is 
located,  i.e.,  cell  0(u)  is  searched  from  vertex  v.  We  denote  the  searcher’s  vertex 
(waypoint)  prior  to  time  period  1  by  uq  ^  V.  This  starting  vertex  could  be  a  desig¬ 
nated  entry  waypoint  into  the  area  of  interest.  We  also  define  V  C  V  to  be  a  set  of 
possible  destination  vertices  for  the  searcher  (e.g.,  V  =  {uq}  if  the  searcher  is  required 
to  return  to  the  starting  vertex  or  V  =  V  if  the  search  can  end  anywhere). 

For  any  k  E  T  and  Vi  E  V,  I  =  0,1,2, k,  such  that  {vi-i,vi)  E  £  for  all 
I  =  l,2,...,k,  let  the  sequence  denote  a  directed  VQ-Vk  subpath.  If  Vk  E  V, 

then  the  directed  VQ-Vk  subpath  is  a  directed  Vo-Vk  path.  When  the  meaning  is  clear 
from  the  context,  we  refer  to  a  directed  vo-Vk  (sub)path  as  a  (sub)path.  In  this 
notation,  the  searcher  flies  from  vq  to  some  Vk  E  V  along  a  directed  Vo-Vk  path.  The 
searcher  occupies  only  one  vertex  v  E  V  each  time  period,  and  stays  at  the  same 
vertex  or  moves  to  another  vertex  in  lF{v)  for  the  next  time  period. 

We  adopt  the  following  target-detection  model.  If  the  target  is  in  cell  c  E  C 
during  time  period  t  E  T  and  the  searcher  is  at  the  same  time  at  waypoint  v'  E  V 
above  cell  c,  i.e.,  (/>(u')  =  c,  then  detection  occurs  with  a  probability.  We  term  this 
probability  as  glimpse  detection  probability  and  refer  to  g{v,v' ,t),  where  u  G  V  is  the 
searcher’s  waypoint  during  time  period  f  — 1.  Hence,  the  glimpse  detection  probability 


11 


during  time  period  t  depends  on  the  previous  and  current  waypoints  for  the  searcher. 
This  is  a  generalization  of  earlier  models  where  the  glimpse  detection  probability 
depends  only  on  v'  and  t  [52,  32]  and  is  one  of  our  contributions.  Our  model  accounts 
for  the  fact  that  moving  from  some  waypoint  to  a  new  waypoint  may  result  in  a  lower 
glimpse  detection  probability  than  if  the  searcher  already  was  loitering  at  the  latter 
waypoint,  i.e.,  g{v',v',t)  >  g{v,v',t)  for  v'  ^  v.  This  effect  occurs  if  refocusing  the 
sensor  and  becoming  familiar  with  a  new  cell  have  a  signihcant  detrimental  effect  on 
the  glimpse  detection  probability.  In  general,  change  of  waypoint,  especially  change 
of  altitude  and  frequent,  irregular  change  of  direction,  may  distract  from  the  search. 
This  generalization  also  allows  us  to  account  indirectly  for  small  transit  times  (much 
less  than  the  length  of  a  time  period)  between  waypoints  without  adopting  a  hne  time 
discretization  with  resulting  high  computational  cost.  In  this  notation,  the  probability 
of  detection  at  waypoint  v'  during  time  period  t,  given  search  at  waypoint  v  during 
time  period  t  —  1,  and  no  prior  detections  becomes  p{(f){v'),t)g{v,v' ,t). 

The  glimpse  detection  probability  may  also  depend  on  the  searcher’s  speed. 
To  account  for  this  situation,  edges  can  be  duplicated  as  in  [10]  to  represent  search 
at  different  speeds  in  cases  where  speed  influences  the  searcher’s  effectiveness  signif¬ 
icantly.  For  the  sake  of  simplicity,  however,  we  do  not  introduce  notation  to  handle 
that  situation. 

Since  p{-,t)  is  the  undetected  target  distribution  at  the  beginning  of  time 
period  t,  it  depends  on  searches  up  to  time  period  t  —  1.  Specihcally,  if  p{-,t)  = 
[p{l,t), . . .  ,p{c',t), . . .  ,p{C,t)]  and  cell  c'  is  searched  from  waypoint  v'  (i.e.,  0(n')  =  c') 
during  time  period  t,  the  undetected  target  distribution  at  the  beginning  of  the  next 
time  period  t  -|-  1  is 

p{-,t  +  l)  =  [p{l,t),..,p{c-l,t),p{c,t){l-g{v,v',t)),p{c+l,t),..,p{C,t)]r,  (II.l) 
where  v  is  the  searcher’s  vertex  during  time  period  t  —  1.  Thus,  the  probability  of 
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detection  along  a  directed  Vo-Vk  path  V  =  denoted  q(V),  is  defined  as 

k 

Qi'P)  =  ^P{<P{vt),t)g{vt-i,vt,t).  (II.2) 

t=i 

Let  X  =  {1,2,...,/}  be  a  set  of  resources  and  fi{v,v\t)  be  the  amount  of 
resource  i  E  I  consumed  by  the  searcher  at  vertex  v'  during  time  period  t,  given  search 
at  vertex  v  during  time  period  t  —  1.  Resources  may  represent  physical  commodities 
such  as  fuel  and  ordnance  that  are  depleted  during  the  search  as  well  as  abstract 
factors  such  as  notions  of  risk  exposure  during  the  search.  In  contrast  to  search  by 
manned  aircraft,  where  significant  risks  are  usually  avoided,  planners  accept  higher 
risks  for  UAV  search  missions  and  would  like  to  balance  risk  with  other  factors  during 
the  planning  process.  We  discuss  the  calculation  of  risk  in  Section  D.  The  total 
“consumption”  of  resource  i  eT  along  the  directed  path  V  =  (wzjjLo 

t).  (II.3) 

t=i 

The  searcher  cannot  consume  more  than  fj  of  resource  i  E  I  along  a  path.  Hence, 
the  resource-constrained  single  searcher  problem  (RSSP)  is  the  problem  to  find  a 
directed  v^-Vk  path  V  =  (vzjjLo)  with  Vk  E  V,  k  E  T,  that  maximizes  q(V)  subject  to 
the  constraints 

TiiV)  <fi,i  El.  (II.4) 

We  contribute  that  the  single  searcher  problem  (SSP)  is  generalized  to  the  case 
with  multiple  resource  constraints.  We  refer  to  constraints  (II. 4)  as  side  constraints. 
In  Section  D,  we  examine  a  case  study  with  two  time-invariant  resources:  risk  and  fuel. 
The  next  section  deals  with  the  “unconstrained”  problem  with  no  side  constraints 
(referred  as  SSP),  while  Sections  C  and  D  address  the  full  RSSP. 

B.  BRANCH-AND-BOUND  ALGORITHM  FOR  SSP 

In  this  section,  we  consider  the  single  search  problem  (SSP)  that  maximizes 
(II. 2)  without  the  side  constraints  (II. 4).  Several  branch- and-bound  algorithms  for 
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the  solution  of  SSP  have  been  developed  [46,  50,  18,  34,  13,  52,  32],  These  algorithms 
implicitly  enumerate  all  searcher  paths  that  cannot  be  proven,  by  means  of  a  bound,  to 
be  nonimproving.  The  next  subsection  presents  the  algorithm  in  [32],  which  appears 
to  be  the  fastest  in  the  literature  for  SSP  under  a  broad  range  of  conditions;  see  [32] 
for  an  empirical  stndy.  The  snbsequent  snbsection  presents  fonr  modihcations  of  that 
algorithm  which  generalize  and  speed  up  that  algorithm. 

This  section  ignores  side  constraints  and,  without  loss  of  generality,  assnmes 
that  an  optimal  path  consists  of  T  +  1  vertices.  To  simplify  the  notation,  we  also 
assume  in  this  section  that  there  is  no  end-point  restriction,  i.e.,  V  =  V.  We  hnd 
identical  assumptions  in  [52,  32].  Initially,  we  assnme  that  the  glimpse  detection 
probability  g{v,v',t)  is  independent  of  v  and  write  g{v',t),  but  relax  that  assumption 
later  in  this  section. 

1.  Existing  Algorithm 

For  completeness  and  ease  of  reference  later,  we  ontline  the  algorithm  in  Lau 
et  ah  [32].  Given  a  subpath  {viYiZq,  t  E  T,  let  /C(t)  be  the  set  of  triplets  of  the  form 
q{vt,t))  representing  extensions  of  yet  to  be  explored.  The  first  element 

Vt  refers  to  the  next  vertex  to  visit,  the  second  element  t  is  the  time  period^  to  visit 
the  vertex  Vt,  and  the  third  element  q{vt,t)  is  an  npper  bonnd  on  the  probability 
of  detection  along  any  path  that  starts  with  the  subpath  {vi}Yo-  The  upper  bound 
qivtit)  consists  of  three  parts.  Let  dt{vt,t)  be  an  upper  bound  on  the  probability  of 
detection  during  time  periods  t-|-l,t-|-2,...,T,  given  that  the  searcher  starts  at  Vt  at 
time  t,  and  no  detection  occnrs  along  the  subpath  {vi}Yo-  The  two  other  parts  are 
the  probability  of  detection  on  the  subpath  and  the  probability  of  detection 

dnring  t.  Hence, 

q{vt,t)  =  q{{vi}\zl)  +  p{(t){vt),t)g{vt,t)  +  dt{vt,t).  (II.5) 

^This  information  is  currently  redundant  but  the  notation  is  convenient  in  later  generalizations. 
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We  also  let  q  denote  the  largest  detection  probability  found  so  far  among  all  the 
examined  paths.  In  this  notation,  the  algorithm  in  [32]  takes  the  following  form. 

Algorithm  1  (Solves  SSP). 

Input:  A  time-expanded  network  {J\f,A)  (dehned  later  as  Figure  3)  with 
nodes  n  E  Af  and  arcs  (n, n')  G  A  where  n  =  {v,t  —  1),  n'  = 

Arc  lengths  g{v',  t)  in  (A/”,  A),  the  initial  target  distribution  p(-,  1),  Markov 
transition  matrix  F,  and  upper  bound  q. 

Output:  An  optimal  search  path  V*  =  {vi}f^Q  and  its  value  q*. 

Step  0.  Set  t  =  0,  IC{t)  =  {(uq,  0,  cxd)},  and  q  =  0. 

Step  1.  If  )C(t)  is  empty,  replace  t  hj  t  —  1.  Else,  go  to  Step  3. 

Step  2.  If  t  =  0,  stop:  the  last  saved  path  is  optimal  and  q  is  its  probability 
of  detection.  Else,  go  to  Step  1. 

Step  3.  Remove  from  IC(t)  the  triplet  {vt,t,q{vt,t))  with  the  largest  bound 
Q{vut). 


Step  4.  If  q{vt,t)  <  q,  go  to  Step  1.  (Current  subpath  is  fathomed.) 

Step  5.  If  t  <  T,  then  for  each  vertex  v  G  A(vt),  calculate  a  bound  dt+i(v,t  + 

1)  as  well  as  q(v,  f  -|-  1)  using  equation  (II. 5),  and  add  (v,  t  +  1,  q{v,  t  +  1)) 
to  Kit  +  1).  Replace  thy  t  +  1  and  go  to  Step  3.  Else,  let  q  =  q{vt,  t)  and 
save  the  incumbent  path  {vi}JLq,  and  go  to  Step  1. 

Clearly,  a  tight  bound  dt{vt,t)  will  reduce  the  number  of  branching  attempts 
in  Step  4  of  Algorithm  1.  As  examined  in  [51,  52,  32],  there  is  a  fundamental  trade-off 
between  the  effort  needed  to  compute  a  bound  and  its  tightness.  From  these  studies, 
it  appears  that  the  bounding  technique  in  [32]  compares  favorably  in  most  situations. 
We  describe  that  bounding  technique  in  the  rest  of  this  subsection. 

Consider  a  subpath  {viYi^q,  t  E  T,  and  let  Pg{-,t)  be  the  undetected  target 
distribution  after  search  along  {vi}Yo,  i-®-) 

Pg{-,t)=  (II.6) 

[p(l,t),..,p(0(nt)  -  l,t),p{(f){vt),t){l  -  g{vt,t)),p{(f){vt)  +  l,t),..,p{C,t)]. 
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We  use  subscript  g  to  indicate  that  Pg{-,t)  is  obtained  from  p{-,t)  by  applying  the 
glimpse  detection  probability  corresponding  to  the  last  vertex  in  the  “current”  sub¬ 
path  {viYi^Q.  For  any  integer  s  >  t,  s  E  T,  we  also  dehne 

pr{;s;t)=Pg{;t)r-\  (II.7) 

As  seen,  pr{c,s-,t)  is  the  probability  that  the  target  is  in  cell  c  at  time  period  s  >  t 
and  there  was  no  detection  during  search  along  the  subpath  {vi}Yo-  Hence,  target 
distribution  pr{-,s]t)  at  time  period  s  >  t  ignores  the  effect  of  search  after  time 
period  t.  If  the  subpath  is  {uq},  i.e.,  t  =  0,  we  dehne  for  notational  convenience 

pr(-,5;0)=p(-,l)r-\  (II.8) 

for  any  s  >  0,  s  G  T.  Moreover,  we  dehne  pr(c,  t;t)  =  0  for  all  c  G  C  and  t  =  0, 1, ...,  T. 

Now,  we  construct  a  time-expanded  graph  from  the  network  (V,£^)  as  follows 
(see  Figure  3).  Each  vertex  u  G  V  is  duplicated  T  times  to  dehne  the  nodes  {u,<s), 
s  E  T.  Let  Af  be  the  set  of  all  such  nodes  as  well  as  the  nodes  uq  =  {vq,  0)  and 
n  =  {v,T+l)  representing  the  searcher’s  prior  position  and  hnal  position,  respectively. 
Here,  v  is  an  artihcial  terminal  vertex.  Two  nodes  n  =  {v,s  —  1)  and  n'  =  {v',s), 
v,v'  eV  and  s  =  2,  3, ...,  T,  are  connected  with  an  arc  (n,  n')  if  and  only  if  (u,  v')  E  S. 
Moreover,  the  node  uq  =  (uo,0)  is  connected  with  an  arc  to  a  node  n'  =  (u',  1), 
v'  E  V,  if  and  only  if  (uo,F)  G  S]  and  every  node  n  =  {v,T),  u  G  V  is  connected 
with  an  arc  to  h.  Let  A  be  the  set  of  all  arcs.  For  any  integer  k  <  T  +  1  and  nodes 
rii  =  {vi,  1)  E  Af,  /  =  0, 1, ...,  k,  such  that  {rii-i,  ni)  E  A  for  all  I  =  1,2,  ...,k,  we  let  the 
sequence  denote  a  subpath  in  the  time-expanded  graph  {Af,A). 

For  some  t  E  {0, 1,  ...,T  —  1},  suppose  that  a  subpath  {vi}Yo  original 

graph  {y,£)  is  given.  Then,  we  endow  each  arc  {n,n')  =  {{v,s),{v',s  +  1))  G  A, 
s  =  t,t  +  l,...,T  —  1,  in  the  time-expanded  graph  {Af,  A)  with  a  “reward” 

Cn,n'  =  [PT{(t>{v'),  s  +  l-,t)  -  s;  t)g{v,  s)F(u,  v')]g{v',  s  -f  1),  (11.9) 

where  F(u,u')  is  the  (j){v)-(j){v')  element  of  the  Markov  transition  matrix  F.  We  set 
Cn,h  =  0  for  all  (n,  h)  E  A.  We  observe  that,  multiplied  out,  the  hrst  term  pv{(t>{v'),  s-f 
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r+1 


0  1  ■■■  s-l  s  5  +  1  ■■■  T 


Figure  3.  A  time-expanded  graph  from  the  network  in  Figure  1  and  one  altitude.  The 
searcher’s  initial  position  is  Uq  =  (uq,  0)  =  ((1, 1),  0)  and  hnal  position  is  h  =  (h,  T-l-1). 

l]t)g{v',s  -|-  1)  in  (II. 9)  is  effectively  the  expected  number  of  detections  during  time 
period  s  -|-  1,  which  gives  rise  to  the  so-called  mean  bound  [34],  and  the  second  term 
in  (II. 9)  improves  the  bound  by  accounting  for  the  effect  of  search  during  time  period 
s  [32].  We  refer  to  (A/*,  .4,)  with  arc  rewards  given  by  (II. 9)  as  the  time-expanded 
network. 

In  Lau  et  al.  [32],  it  is  shown  that  given  the  subpath  the  optimal  value 

of  the  longest-path  problem  in  the  time-expanded  network  from  node  (ut,t)  to  node 
('0,T  -|-  1),  using  the  rewards  in  (II. 9)  as  “arc  length,”  provides  an  upper  bound  for 
Algorithm  1.  Specihcally,  this  optimal  value  is  an  upper  bound  on  the  probability  of 
detection  during  time  periods  t-|-l,t-|-2,...,T,  given  that  the  searcher  starts  at  Vt  at 
time  t,  and  no  prior  detections  occurred  along  the  subpath  We  denote  this 

bound  by  dt{vt,  t)  and  refer  to  it  as  the  dynamic  bound,  as  it  needs  to  be  recomputed 
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every  time  the  current  subpath  is  extended. 

Since  the  time-expanded  network  is  acyclic,  the  longest-path  problem  can  be 
solved  in  polynomial  calculation  time  with  a  standard  shortest-path  algorithm  ([1], 
pages  77-79).  We  observe  that  to  compute  dt{vt,t)  given  the  subpath  it  is 

only  necessary  to  generate  the  part  of  the  graph  {J\f,A),  and  the  corresponding  arc 
rewards,  “after”  time  t  and  within  reach  from  node  {vt,t),  since  the  longest  path 
starts  at  node  {vt,t)-  Hence,  in  Step  5  of  Algorithm  1,  it  suffices  to  generate  the 
time-expanded  network  “after”  time  t. 

2.  Algorithmic  Modifications  for  SSP 

Algorithm  1  requires  bound  calculation  at  each  branching  which  is  a  principal 
bottleneck  for  that  algorithm.  Thus  we  present  new  bounding  techniques  to  enhance 
Algorithm  1.  This  subsection  proposes  and  examines  four  modifications  of  Algorithm 
1.  The  hrst  modification  extends  the  bound  in  [32]  to  the  case  in  which  glimpse 
detection  probability  depends  on  the  previous  vertex,  i.e.,  g{v,v',t).  The  following 
three  modihcations  are  aimed  to  simplify  and  speed  up  Algorithm  1.  The  second 
modihcation  weakens  the  bound  but  greatly  simplihes  its  bound  calculation.  The 
third  modihcation  improves  the  weaker  bound  at  little  computational  expense.  The 
fourth  modihcation  takes  advantage  of  a  special,  but  frequently  occurring,  initial  tar¬ 
get  distribution. 


a.  Bound  for  Edge-Dependent  Glimpse  Detection  Proba¬ 
bility 

Previous  studies  (see,  e.g.,[52,  32])  assume  that  the  glimpse  detection 
probability  depends  only  on  the  searcher’s  current  waypoint  (vertex)  and  on  time. 
As  we  argue  in  Section  A,  this  is  somewhat  restrictive.  Fortunately,  we  can  easily 
extend  the  previous  studies  to  the  case  where  the  glimpse  detection  probability  also 
depends  on  the  searcher’s  previous  waypoint.  The  only  modihcation  that  is  required 
is  to  redehne  the  arc  reward  Cny  in  (II. 9).  However,  a  straightforward  replacement  of 
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g{v,  s)  by  g{y" ,  v,  s)  in  (II. 9),  where  v”  is  the  vertex  prior  to  v,  would  ruin  the  longest- 
path  structure  of  the  bound-calculation  problem;  Cn,n'  would  no  longer  depend  only  on 
the  head  and  tail  of  the  arc  (n,  n').  Hence,  it  is  necessary  to  use  the  smallest  glimpse 
detection  probability  min„//g7^(„)  g{v'',  v,  s)  to  eliminate  the  dependence  on  the  vertex 
prior  to  v,  where  7^(w)  C  V  is  the  reverse  star  of  v,  i.e.,  7^(v)  =  {v"  G  V  |  (v",  v)  G  £}. 
Consequently,  we  now  endow  each  arc  (n,  n')  =  {{v,  s),  {v',  s-|-l))  G  A  with  the  reward 


Cnri/  — 


Pr{<l){v),s  +  l;t)  -pr(0(v),s;f)  min  g{v',v,s)  T{v,v) 

\v"e:Tl{v)  / 


g{v,v',s  +  l). 

(11,10) 

With  this  reward,  the  bound  calculation  remains  a  longest-path  problem  in  an  acyclic 
graph  and  it  can  be  shown  using  the  same  arguments  as  in  [32]  which  prove  that  the 
dynamic  bound  is  valid. 


b.  Static  Bound 

Algorithm  1  requires  one  longest-path  calculation  in  the  time-expanded 
network  for  each  vertex  in  the  forward  star  of  the  current  vertex  to  compute  the  re¬ 
quired  bounds  dt{vt,t)  (see  Step  5).  The  approach  of  Algorithm  1  with  dynamic 
bound  follows  the  traditional  approach  of  branch-and-bound  algorithms  where  the 
bound  is  reoptimized  before  each  branching.  In  the  present  case,  the  reoptimization 
corresponds  to  the  longest-path  calculation  and  requires  computing  the  arc  rewards 
Cn^n',  see  (II.  10).  These  calculation  can  be  time  consuming,  as  they  typically  involve 
a  moderately  large  Markov  transition  matrix  T  and  associated  matrix  multiplication. 
Our  numerical  example  considers  that  T  is  of  size  225  by  225.  We  propose  to  use  a 
static  bound  instead  of  the  dynamic  bound  proposed  in  [32]  and  described  in  Subsec¬ 
tion  B.l.  As  shown  below,  all  the  necessary  static  bounds  are  computed  prior  to  any 
branching  and  are  not  recomputed  later. 

The  dynamic  bound  dt{vt,t)  (see  Subsection  B.l)  uses  information 
about  search  along  a  current  subpath  {viYi^q.  However,  the  bound  remains  valid 
if  the  effect  of  search  along  the  current  subpath  is  ignored.  This  follows  from  the 
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same  arguments  as  in  the  proof  of  the  validity  of  dt{vt,t),  see  [32],  We  denote  the 
new  bound  do{vt,t),  where  the  subscript  0  indicates  that  the  trivial  subpath  {uq} 
is  used  in  (II. 10)  with  t  =  0  instead  of  the  subpath  which  dt{vt,t)  utilizes. 

Hence,  the  only  difference  between  dt{vt,  t)  and  d^ivt,  t)  is  that  pr(‘,  S  i)  is  replaced  by 
Pr(',s0)  ill  (IhlO).  However,  this  difference  makes  the  bound  do{vt,t)  independent 
of  the  current  subpath  used  to  reach  the  vertex  u*.  Hence,  do{v,t)  can  be  computed 
in  advance  for  all  nodes  {v,t)  G  A/”,  and  dynamical  computation  of  bounds  is  not 
required.  Consequently,  the  arc  rewards  (11.10)  and  bounds  are  computed  only  once. 
We  refer  to  do{v,t)  as  the  static  bound.  We  observe  that  it  is  not  necessary  to  carry 
out  a  longest-path  calculation  from  each  node  {v,  t)  G  A/”  to  (h,  T  +  1)  to  obtain 
do{v,t).  It  is  more  efficient  to  carry  out  the  longest-path  calculations  backward  from 
node  (h,  T  -|-  1)  to  all  nodes.  This  calculation  simply  amounts  to  applying  a  shortest- 
path  algorithm  once  to  the  time-expanded  network  with  arc  lengths  equal  to  the 
negative  rewards. 

In  Step  5  of  Algorithm  1,  we  now  simply  use  do{vt,  t)  instead  of  dt{vt,  t). 
Thus,  the  modihed  algorithm  does  not  require  any  longest-path  calculation  in  Step  5. 
All  bound  calculations  are  done  prior  to  Step  0.  Clearly,  the  modihed  approach  results 
in  a  weaker  bound  and  the  need  for  more  branching  attempts.  However,  the  additional 
branching  attempts  may  be  compensated  by  shorter  per-iteration  computing  times. 

In  order  to  examine  the  effect  of  the  static  bound  do{vt,t),  we  examine 
the  same  numerical  example  as  in  [32]:  An  area  of  interest  consists  of  11  by  11  cells. 
The  searcher  operates  only  at  one  altitude  and  its  moves  are  restricted  to  vertically 
and  horizontally  adjacent  cells,  excluding  diagonal  moves.  The  target  remains  in  the 
current  cell  with  a  probability  p  or  moves  to  one  of  the  vertically  or  horizontally 
adjacent  cells  with  probability  1  —  p.  The  different  moves  are  equally  likely.  The 
searcher  departs  cell  1  (uq  =  1),  where  the  cells  are  numbered  from  left  to  right  and 
from  top  to  bottom.  Hence,  cell  1  is  the  upper-left-corner  cell.  The  target  starts  at 
the  center  cell,  i.e.,  at  cell  61.  The  time  horizon  T  =  17. 
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We  have  implemented  Algorithm  1  with  static  bound  using  Microsoft 
Visual  C++  6.0  on  a  desktop  computer  with  a  3.4  GHz  Intel  Pentium  IV  processor, 
1.0  gigabytes  of  RAM,  and  the  Microsoft  Windows  XP  operating  system.  Table  1 
shows,  for  a  range  of  constant  glimpse  detection  probabilities  g{v,v',t)  and  “stay 
probabilities”  p,  the  run  times  (in  seconds)  and  numbers  of  bounding  attempts  of 
Algorithm  1,  as  well  as  those  reported  in  [32].  We  especially  consider  the  case  of  high 
glimpse  detection  probability  {g{v,v',t)  =  0.99),  in  which  both  static  and  dynamic 
bounds  tend  to  be  relatively  weak.  Column  3  of  Table  1  presents  the  resulting  run 
times  for  Algorithm  1  with  static  bound.  The  corresponding  run  times  reported  from 
[32]  are  found  in  column  5.  (The  case  g{v,v',t)  =  0.99  is  not  considered  in  [32].) 
Those  reported  numbers  are  achieved  on  a  2.6  GHz  Opteron  152  processor  computer 
using  Matlab.  Hence,  a  direct  comparison  between  the  run  times  in  columns  3  and 
5  is  difficult.  However,  we  hnd  reports  of  run  times  for  other  problem  instances 
using  a  C++  implementation  in  [32].  A  brief  comment  in  [32]  based  on  a  single 
problem  instance  indicates  that  the  Matlab  implementation  is  15.6  times  slower  than 
the  C++  implementation.  For  the  sake  of  comparison,  we  scale  down  the  run  times 
in  column  5  of  Table  1  with  1/15.6  to  approximately  account  for  the  slower  Matlab 
implementation.  We  also  scale  down  the  run  times  of  column  5  by  a  factor  of  2. 6/3. 4 
to  account  for  the  slower  computer  used  in  [32].  The  resulting  scaled  down  run  times 
are  presented  in  column  6  of  Table  1.  We  observe  that  the  scaled  down  run  times  for 
Algorithm  1  with  dynamic  bound  are  somewhat  faster  than  the  corresponding  run 
times  with  static  bound.  However,  the  static  bound,  which  is  weaker  than  the  dynamic 
bound,  remains  fairly  competitive,  especially  for  more  difficult  problem  instances. 

We  compare  the  strengths  of  the  static  and  dynamic  bounds  by  counting 
the  number  of  branching  attempts  required  in  Algorithm  1.  Columns  4  and  7  of 
Table  1  give  the  numbers  of  branching  attempts  for  the  static  and  dynamic  bounds, 
respectively.  The  number  of  branching  attempts  for  the  static  bound  is,  on  average, 
37  times  larger  than  in  the  case  of  the  dynamic  bound.  We  observe  that  the  greater 
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Static  Bound 

Dynamic  Bound  [32] 

Dynamic  Bound 

Time 

Branching 

Time 

Scaled 

Branching 

Time 

Branching 

g{v,v',t) 

P 

(sec.) 

attempts 

(sec.) 

(sec.) 

attempts 

(sec.) 

attempts 

0.3 

3.59 

2,301,182 

14.56 

0.71 

58,314 

2.09 

58,314 

0.3 

0.6 

2.58 

1,635,517 

14.53 

0.71 

52,369 

2.11 

52,369 

0.9 

11.47 

7,338,492 

157.30 

7.71 

380,889 

20.75 

380,889 

0.3 

5.38 

3,424,282 

19.07 

0.94 

49,774 

2.59 

49,774 

0.6 

0.6 

2.58 

1,620,402 

23.76 

1.16 

47,454 

3.03 

47,454 

0.9 

59.26 

38,186,809 

730.96 

35.84 

2,185,066 

103.08 

2,185,066 

0.3 

5.78 

3,675,197 

21.55 

1.06 

45,019 

2.78 

45,019 

0.9 

0.6 

2.89 

1,831,875 

30.58 

1.50 

59,527 

3.91 

59,527 

0.9 

176.26 

113,646,357 

2902.27 

142.30 

11,299,259 

431.38 

11,299,259 

0.3 

6.09 

3,865,002 

N/A 

N/A 

N/A 

2.91 

46,339 

0.99 

0.6 

3.00 

1,896,960 

N/A 

N/A 

N/A 

3.91 

58,752 

0.9 

192.06 

123,822,672 

N/A 

N/A 

N/A 

558.63 

16,685,969 

Table  1.  Run  times  and  number  of  branching  attempts  (counted  in  Step  4)  for 
Algorithm  1  with  static  and  dynamic  bounds  on  11  by  11  cell  search  problem  with 
time  horizon  T  =  17.  Columns  labeled  “Dynamic  Bound  [32]”  correspond  to  original 
and  speed- adjusted  results  from  [32], 


number  of  bounding  attempts  is  partially  compensated  for  by  avoiding  dynamical 
reoptimization  of  the  bound. 

For  a  direct  comparison  between  the  static  and  dynamic  bounds,  we 
implement  Algorithm  1  with  dynamic  bound  in  Microsoft  Visual  C-I--I-  6.0.  We  made 
a  signihcant  effort  to  ensure  that  the  implementation  is  efficient,  including  efficient 
handling  of  sparse  matrices.  Column  8  and  9  of  Table  1  report  the  run  times  and  the 
number  of  branching  attempts  for  our  implementation  of  Algorithm  1  with  dynamic 
bound,  respectively.  We  observe  that  our  implementation  results  in  identical  numbers 
of  branching  attempts  compared  to  the  implementation  in  [32].  When  comparing 
columns  6  and  8,  we  see  that  our  implementation  of  the  dynamic  bound  results  in 
somewhat  longer  run  times  than  the  scaled  times  from  [32] .  However,  the  longer  times 
in  column  8  compared  to  column  6  can  partially  be  due  to  an  excessively  aggressive 
scaling  of  run  times  going  from  column  5  to  column  6. 

While  implementing  the  dynamic  bound,  we  noted  a  signihcant  chal¬ 
lenge  associated  with  efficient  matrix  multiplication  and  data  handling.  Since  the 
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dynamic  bound  dt{v,t)  is  reoptimized  a  large  number  of  times,  it  is  paramount  to 
carry  out  the  related  calculations  efficiently.  In  comparison,  it  is  rather  trivial  to 
implement  a  static  bound.  It  is  not  critical  to  carry  out  the  longest-path  calculations 
highly  efficiently  in  the  static  case  as  they  are  only  done  once. 

We  also  observe  that  both  static  and  dynamic  bounds  tend  to  be  weaker 
when  the  target  is  nearly  stationary  (e.g.,  p  =  0.9)  and  glimpse  detection  probability 
is  large  (e.g.,  g{v,v',t)  =  0.9  or  0.99).  For  these  problem  instances,  the  implementa¬ 
tion  with  dynamic  bound  is  more  sensitive  to  the  change  in  data.  In  fact,  comparing 
the  instance  {g{v,v',t)  =  0.9  and  p  =  0.9)  to  the  instance  {g{v,v',t)  =  0.99  and 
p  =  0.9),  we  see  that  the  run  time  and  the  number  of  branching  attempts  in  the 
case  of  the  dynamic  bound  become  1.29  and  1.48  times  larger,  respectively.  On  the 
other  hand,  the  static  bound  is  less  sensitive,  and  its  numbers  become  only  1.09  times 
larger.  These  numbers  indicate  that  it  becomes  even  less  worthwhile  to  invest  time 
in  dynamic  bound  calculations  when  the  bounds  are  relatively  weak  anyway. 

c.  Directional  Static  Bound 

As  seen  from  Table  1,  the  static  bound  is  substantially  weaker  than 
the  dynamic  bound.  We  derive  a  stronger  static  bound  motivated  by  the  classical 
approach  to  handling  turn-radius  constraints  in  vehicle-routing  problems  [8] . 

In  the  longest-path  calculations  for  the  static  bound,  the  reward  of  arc 
{{v,s),  {v',s  +  1))  is,  effectively,  the  probability  of  detection  at  vertex  F  during  time 
period  s  -|-  1  and  no  detection  at  vertex  v  during  time  period  s.  Of  course,  this 
overestimates  the  probability  of  detection  at  vertex  v'  during  time  period  s  -|-  1  and 
no  prior  detections,  as  detection  could  occur  prior  to  time  period  s.  We  strengthen 
the  static  bound  if  we  redehne  the  arc  reward  to  be  the  probability  of  detection  at 
vertex  v'  during  time  period  s  -|-  1  and  no  detection  at  vertex  v  during  time  period  s 
and  no  detection  at  the  vertex  visited  during  time  period  s  —  1.  However,  redehning 
the  arc  reward  to  depend  not  only  on  the  arc’s  head  and  tail  nodes,  but  also  on  a 
previous  node,  ruins  the  longest-path  structure  of  the  bound-calculation  problem. 
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A  similar  situation  arises  in  vehicle  routing  problems  for  vehicles  with 
turn-radius  constraints  or  penalties.  The  classical  approach  to  handle  that  situation 
is  to  duplicate  each  node  a  number  of  times  equal  to  the  number  of  nodes  in  the 
node’s  reverse  star.  An  arc  in  the  resulting  “node-expanded”  network  then  carries 
information  about  three  nodes,  not  only  two,  and  a  desirable  network  structure  of 
the  problem  can  be  maintained.  Fortunately,  it  is  practical  to  carry  out  such  a  node¬ 
expansion  approach  in  the  problems  of  interest  in  this  paper  because  the  number  of 
nodes  in  the  reverse  star  is  typically  quite  moderate.  Hence,  we  proceed  along  the 
stated  lines  and  develop  a  node-and-time  expanded  network,  in  which  the  improved 
static  bound  can  be  calculated  by  solving  a  longest-path  problem.  We  refer  to  this 
improved  bound  as  the  directional  static  bound. 

For  any  n'  G  A/”,  let  lZ{n')  d  N  he  the  reverse  star  of  n',  i.e.,  lZ{n')  = 
{n  G  N'\{n,n')  G  A}.  Then,  for  any  n,n'  G  J\f\{h}  such  that  (n, n')  G  A,  we  define 
an  expanded  node  ^  =  {n,n').  We  do  not  expand  the  end  node,  so  we  set  ^  =  n.  Let 
S  be  the  set  of  all  expanded  nodes.  Two  expanded  nodes  G  S  are  connected  by 
an  expanded  arc  if  ■C  =  {n,n')  and  =  {n',n").  Let  the  set  of  all  expanded 

arcs  be  H.  The  node-and-time  expanded  graph  is  illustrated  in  Figure  4. 

We  endow  each  expanded  arc  in  the  node-and-time  expanded  graph 
(S,  H)  with  a  reward  similar  to  (11.10).  To  derive  the  exact  form  of  this  reward,  we 
need  the  following  building  blocks.  For  any  G  V  and  t  G  T,  let  Mt{v,v')  he  a  C 
by  C  identity  matrix  with  the  0(n')-th  diagonal  element  set  equal  to  1  —  g{v,v',t)- 
We  also  let  F(n')  be  the  column  of  the  Markov  transition  matrix  F. 

From  (11.2)  and  the  recursive  application  of  (ILl),  we  see  that  the 
probability  of  detection  along  a  path  {vi}f^Q  is  given  by 

q{{vi}f=o)  =  p{(j){vi),  l)g{vo,  ui,  1)  + 
p(-,  l)Mi{vo,vi)T{v2)g{vi,V2,2)  + 
p{-,  l)Mi{vo,  vi)rM2{vi,  V2)T{v3)g{v2,  Vs,  3)  -h  (11.11) 
p(-,  l)Mi(no,  ni)FM2(ni,  n2)FM3(n2,  n3)F(n4)5((n3,  ^4, 4)  4- 
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Figure  4.  A  node-and-time  expanded  network  from  the  time-expanded  graph  (Figure 
3). 


p(-,  l)Mi{vo,Vi)TM2{Vi,V2)  • . . .  ■  TMT-l{vT-2,VT-l)'r{vT)g{vT-l,VT,T). 

The  expression  (IF  11)  gives  insight  into  a  class  of  bounds  on  the  probability  of  detec¬ 
tion,  including  the  static  bound  do{vt,t).  If  we  replace  Mt{-,  •)  by  the  identity  matrix 
in  (II.  11),  we  hnd  that 

q{{vi}f^o)  <  l)g{vo,  ui,  1)  + 

p{-,l)T{v2)givi,V2,2)  + 

p(-,l)FF(n3)^(n2,n3,3)  +  (11.12) 

p{-,  l)FFF(n4)5((n3,n4,4)  -t- 

p(-,l)F^  ‘^T{vT)g{vT-i,VT,T). 
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In  (11.12),  the  “reward”  received  during  a  time  period  is  simply  the  expected  number 
of  detection  during  that  time  period  and  depends  only  on  the  current  and  previous 
vertices.  Hence,  it  is  possible  to  compute  an  upper  bound  on  the  optimal  probability 
of  detection  by  hnding  a  path  {vi}f^Q  that  maximizes  the  right-hand  side  in  (11.12). 
This  calculation  amounts  to  a  longest-path  problem  and  is,  in  fact,  the  approach 
in  [34].  (Note,  however,  that  [34]  assumes  that  the  glimpse  detection  probability  is 
independent  of  the  previous  vertex.) 

If  we  replace  each  Mt{-,  ■)  by  the  identity  matrix  everywhere  except  the 
last  matrix  of  each  line  in  (II.  11),  we  obtain 

q{{vi}f^o)  <  p{(j){vi),  l)g{vo,Vi,  1) 

p{-,  l)Mi{vo,vi)T{v2)g{vi,V2,2)  + 

p(-,  l)TM2{vi,V2)T{v3)g{v2,V3,3)  +  (11.13) 

p{-,  l)TTM3{v2,V3)T{vi)g{v3,V4,  A)  + 


pi':  1)1^"^  ^Mr-i(nT-2,  VT-i)T{vT)givT-i,  vt,  T). 

Now,  the  reward  received  during  each  time  period  also  depends  on  the  searcher’s 
position  two  time  periods  ago,  and  the  problem  of  hnding  a  path  that  maximizes 
the  right-hand  side  is  no  longer  a  longest-path  problem.  However,  the  bound  remains 
valid  with  the  following  minor  modihcation,  where  the  maximization  of  a  matrix  with 
a  single  element  different  from  zero  or  one  is  simply  the  maximization  of  that  element; 

qi{vi}f^o)  <  p(0(ui),  l)g{vo,Vi,  1)  + 
p(-,l)  (  niax  Mi(n,ni)  )  r(n2)^(ni,  ua,  2)  + 

\v£-R{vi)  J 

p(-,l)r(  max  M2(n,n2) )  ^(^;3)5((^;2,^'3,3) -h  (11.14) 

pi',  l)rr  (  max  M3(n,n3)  )  T{v4)g{v3,V4,4:)  + 

\v£TI{v3)  I 
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p(-,l)r'^^(  max  MT-iiv,VT-i)]T{vT)givT-i,VT,T). 

\v£-RivT-i)  J 

After  this  modification,  we  see  that  the  reward  during  each  time  period  only  depends 
on  the  current  and  previous  vertices.  Hence,  again,  it  is  possible  to  compute  an  upper 
bound  on  the  optimal  probability  of  detection  by  solving  a  longest-path  problem.  In 
fact,  this  is  exactly  the  approach  we  described  in  Subsection  B.2a  and  it  can  be  shown 
that  the  reward  in  the  longest-path  problem  see  (II.  10),  can  be  deduced  from 
(11.14).  Specihcally,  when  the  current  subpath  in  (II. 10)  is  {uq},  we  have  for  arc 
{n,n')  =  {{v,  s),  {v',  s  -I-  1))  G  ^  that 

Cn,n' =  (  max  M4n",n) )  r(n')fl'(n,n',s -h  1).  (11.15) 

\v"£TZ{v)  j 

Using  similar  arguments,  we  dehne  the  directional  static  bound  as  fol¬ 
lows.  Clearly, 

q{{vi}f=o)  <  l)g{vo,  vi,  1)  + 

p{-,  l)Mi{vo,vi)T{v2)g{vi,V2,2)  + 

p(-,l)  (  max  Mi(n,ni)  ]  rM2(ni,n2)r(n3)5f(n2,n3,3) -h(II.16) 

\ven{vi)  J 

p(-,l)r(  max  M2{v,V2)]TM2,{v2,V:^)T{vi)g(vz,Vi,A)  + 

\v&'R.{v2)  j 


p(-,l)r^  ^  max  MT-2{v,VT-2)]'^MT-l{vT-2,VT-l)T{vT)g{vT-l,VT,T). 
\v&n{vT-2)  J 

Hence,  we  can  compute  an  upper  bound  on  the  optimal  probability  of  detection  by 
hnding  a  path  maximizes  the  right-hand  side  of  (11.16).  This  calculation 

amounts  to  a  longest-path  problem  in  the  node-and-time  expanded  graph  (S,r2). 
The  arc  reward  in  this  longest-path  problem  is  deduced  from  (11.16).  Specihcally,  an 
expanded  arc  (^,^0  =  ^  with  n"  =  {v'\s  —  1),  n  =  (n,s),  and 

n'  =  (n',  s  -f  1),  is  endowed  with  the  reward 


=  p(-,  i)r 


s-2 


max 


M,_i(U",U')  VM,{v'\v)V{v')g{v,v\s  + 1).  (11.17) 


27 


We  refer  to  the  node-and-time  expanded  graph  (S,  ^2)  with  the  arc  rewards  from 
(11.17)  as  the  node-and-time  expanded  network.  Since  the  node-and-time  expanded 
graph  is  acyclic,  longest-path  problems  are  solvable  by  standard  shortest-path  algo¬ 
rithms. 

In  view  of  the  above  discussion,  we  obtain  the  following  result. 

Proposition  1  For  any  v'  eV  and  t  E  T,  let 

(i)  do{v',t)  be  the  value  of  the  longest-path  from  node  {v',t)  to  node  n  in  the 
time-expanded  graph  {Af,A)  with  are  rewards  given  by  (11.15),  and 

(a)  6o{v,v',t)  be  the  value  of  the  longest-path  from  expanded  node  ((n,t  — 1),  {v' t)) 
to  expanded  node  f  in  the  node-and-time  expanded  graph  (S,  hi)  with  are  re¬ 
wards  given  by  (11.17). 

Then,  both  do{v',t)  and  So{v,v',t)  are  upper  bounds  on  the  probability  of  detection 
during  time  periods  t  -\-  l,t  -\-  2,  ...,T  for  any  path  {vi}JLq  with  Vt-i  =  v  and  Vt  =  v' . 
Moreover,  So(v,v',t)  <  do(v',t). 

We  refer  to  So(v,v',t)  as  the  directional  static  bound  and  see  from 
Proposition  1  that  it  is  at  least  as  strong  as  the  static  bound.  We  demonstrate  in  an 
empirical  study  below  that  it  may  be  substantially  stronger. 

Clearly,  building  the  node-and-time  expanded  graph  (S,r2),  computing 
the  associated  rewards,  and  calculating  the  longest-paths  take  some  computing  time. 
However,  the  process  is  only  carried  out  once  before  the  start  of  Algorithm  1  and 
the  computed  bounds  are  stored  for  later  use.  Hence,  the  time  for  computing  the 
directional  static  bounds  remains  small  compared  to  the  overall  run  time  of  Algorithm 
1.  We  conjecture  that  the  use  of  the  node-and-time  expanded  graph  (2,12)  with 
dynamic  reoptimization  of  bounds  will  not  be  efficient  due  to  the  small  but  signihcant 
effort  required  to  build  the  node-and-time  expanded  graph,  compute  associated  arc 
rewards,  and  carry  out  the  longest-path  calculations.  We  observe  that  this  conjecture 
appears  to  be  aligned  with  [32] ,  which  alludes  to  the  inefficiency  of  a  dynamic  bound 
based  on  more  than  one-time-step  look-behind. 


In  Table  2,  we  report  computational  results  for  Algorithm  1  with  the 
directional  static  bound  applied  to  the  same  problem  instances  as  in  Table  1.  Columns 
3  and  4  give  run  times  and  number  of  branching  attempts,  respectively,  for  Algorithm 
1  using  the  directional  static  bound.  We  observe  that,  on  average,  the  number  of 
branching  attempts  has  been  reduced  to  58.1%  by  using  the  directional  static  bound 
compared  with  the  static  bound.  (Compare  column  4  in  Table  1  with  column  4  in 
Table  2.)  Similarly,  the  run  times  have  been  reduced  to  58.2%  by  using  the  directional 
static  bound.  Since  the  reduction  in  branching  attempts  and  run  time  are  essentially 
identical,  we  conclude  that  the  time  for  computing  the  directional  static  bound  is 
small  compared  to  the  overall  run  time. 


g{v,v',t) 

P 

Algo.  1: 

Time 

(sec.) 

D-Static 

Branching 

attempts 

Algo.  2: 

Time 

(sec.) 

Static  &  Red. 
Branching 
attempts 

Algo.  2: 

Time 

(sec.) 

D-Static  &  Red. 
Branching 
attempts 

0.3 

2.59 

1,642,619 

0.22 

129,990 

0.19 

91,198 

0.3 

0.6 

1.80 

1,125,929 

0.17 

94,307 

0.16 

65,485 

0.9 

6.30 

4,032,951 

0.58 

354,677 

0.41 

223,435 

0.3 

3.50 

2,245,784 

0.31 

194,425 

0.27 

128,526 

0.6 

0.6 

1.45 

902,409 

0.17 

92,997 

0.14 

57,015 

0.9 

29.67 

19,321,387 

2.17 

1,442,612 

1.37 

844,747 

0.3 

3.59 

2,282,989 

0.34 

209,160 

0.28 

138,598 

0.9 

0.6 

1.66 

1,037,829 

0.17 

101,216 

0.17 

61,005 

0.9 

80.50 

52,527,302 

4.95 

3,323,101 

2.89 

1,837,647 

0.3 

3.77 

2,396,764 

0.36 

219,217 

0.28 

137,063 

0.99 

0.6 

1.72 

1,080,459 

0.19 

102,763 

0.17 

62,890 

0.9 

87.99 

57,427,410 

5.39 

3,592,696 

3.09 

1,974,871 

Table  2.  Run  time  and  number  of  branching  attempts  (counted  in  Step  4)  for  Al¬ 
gorithm  1  on  problem  instances  of  Table  1  using  directional  static  bound  (D-Static) 
and  for  Algorithm  2  using  static  bound  and  network  reduction  (Static  &  Red.)  and 
directional  static  bound  and  network  reduction  (D-Static  &  Red.). 


d.  Network  Reduction 

In  some  practical  situations,  the  searcher’s  and  the  target’s  initial  po¬ 
sitions  are  relatively  far  from  each  other.  Hence,  for  a  number  of  time  periods  the 
searcher  will  only  examine  cells  where  the  target  is  guaranteed  not  to  be  located.  In 
these  initial  moves,  the  searcher  is  simply  positioning  itself  for  the  later  search.  This 
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Time  period  1 


Time  period  2 


Figure  5.  An  area  of  interest  composed  of  5  by  5  cells.  For  each  time  period,  the 
unshaded,  lightly-shaded,  heavily-shaded,  and  completely-shaded  cells  describe  the 
regions  where  neither  searcher  nor  target  stay,  only  target  possibly  stays,  only  searcher 
possibly  stays,  and  target  and  searcher  possibly  stay,  respectively. 

is  the  situation  in  the  numerical  examples  of  [32].  In  this  subsection,  we  derive  a 
network  reduction  technique  that  utilizes  this  situation. 

Consider  the  time-expanded  graph  Suppose  that  the  searcher 

flies,  for  the  hrst  s  time  periods,  along  a  subpath  {nt}f=Q,  with  Ut  =  such 

that  p{(l){vt),t)  =  0  for  all  t  <  s,  and  p(0(t's),<s)  >  0.  Then,  the  searcher  cannot 
detect  the  target  prior  to  time  period  s.  We  refer  to  the  last  node  of  the  subpath 
{nt}f^Q  as  a  node-of-£rst-contact.  It  is  typically  straightforward  to  determine  a  set  of 
nodes-of-first-contact  by  applying  a  standard  search  algorithm  ([1],  pages  73-77)  on 
the  time-expanded  network.  In  Figure  5,  we  illustrate  a  situation  where  the  searcher 
can  at  the  earliest  detect  the  target  during  time  period  4,  and  also  note  that  the  time 
period  of  all  nodes-of-£rst-contact  {vg,  s)  E  J\f  is  s  >  4. 
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We  observe  that  there  might  be  several  different  subpaths  rit  = 

{vtit)  to  reach  a  node-of-£rst-contact  {vs,s)  G  J\f.  However,  the  subpath  used  to 
reach  {vs,s)  is  of  no  or  little  importance.  In  fact,  if  g{vs-i,Vs,  s)  does  not  vary 
with  the  choice  of  n^-i,  all  subpaths  to  {vs,s)  have  probability  of  detection  equal 
to  p{(l){vs),  s)g{vs-i,Vs,  s)]  with  p{(l>{vs),s)  =  p(-,  l)r^“^r(ns).  Thus,  it  is  enough  to 
find  a  subpath  to  each  node-of- first-contact  {vg,  s)  using  a  standard  search  algorithm, 
connect  initial  node  {vq,  0)  to  {vg,  s)  directly  with  a  “jump”  arc  representing  the 
move  to  {vs,s),  and  ignore  time  periods  l,2,...,s  —  1  during  the  branch-and-bound 
algorithm.  Clearly,  this  procedure  may  reduce  the  amount  of  branching  attempts 
significantly.  We  can  generalize  this  to  the  case  with  g{vg-i,  Vg,  s)  varying  with  respect 
to  ns_i,  as  discussed  in  Subsection  C.2c. 

If  a  lower  bound  q  on  the  optimal  probability  of  detection  is  available  in 
advance  of  the  network  reduction  procedure  described  above,  it  may  not  be  necessary 
to  consider  all  nodes-of-first-contact.  Fortunately  in  the  single  searcher  problems 
(SSP),  a  lower  bound  is  trivially  obtained  by  computing  the  probability  of  detection 
along  the  path  corresponding  to  the  static  bound  (io('i^o,0).  Furthermore,  for  each 
node-of-hrst-contact  {vg,s),  an  upper  bound  on  the  probability  of  detection  after 
time  period  s  with  no  prior  detections  (e.g.,  the  static  bound  do{vg,s))  is  available 
in  advance  of  both  branch-and-bound  and  network-reduction  procedures.  The  upper 
and  lower  bounds  can  be  used  to  eliminate  some  nodes-of-first-contact  that  cannot 
lie  on  an  optimal  path.  Specifically,  if  p{(l){vs),  s)g{vg_i,Vg,  s)  -|-  do{vg,s)  <  q,  we  can 
eliminate  {vg,  s)  and  all  incoming  and  outgoing  arcs.  Thus,  the  amount  of  branching 
attempts  may  be  reduced  further. 

In  order  to  take  advantage  of  this  reduced  network  in  the  branch-and- 
bound  algorithm,  we  construct  a  second  algorithm  (Algorithm  2)  that  generalizes 
the  branch-and-bound  mechanism  of  Algorithm  1.  Before  we  describe  the  algorithm, 
we  need  to  clarify  some  notation.  After  the  algorithm  is  presented,  we  describe  the 
network  reduction  procedure  in  detail. 
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Given  a  subpath  {nj}|hQ,  ni  =  k  E  T,  let  /C(i),  as  before,  be 

the  set  of  triplets  {vt,t,q{vt,t))  representing  extensions  of  yet  to  be  explored. 

Now,  the  index  i  of  /C(i)  refers  to  the  depth  from  node  {vo,  0)  to  the  next  node  {vt,t) 
in  the  branch-and-bound  tree.  Since  (no^O)  is  directly  connected  to  each  node-of- 
hrst-contact  {vs,s),  the  corresponding  depth  in  the  branch-and-bound  tree  is  1.  For 
reporting  purposes,  a  subpath  is  stored  for  each  node-of-£rst-contact  {vs,s). 

Algorithm  2  (Solves  SSP). 

Input:  A  time-expanded  network  with  arc  length  Cn^  or  a  node-and- 

time  expanded  network  (S,  G)  with  arc  length  The  initial  target 

distribution  p{-,  1),  Markov  transition  matrix  F,  and  upper  bound  q. 

Output:  An  optimal  search  path  V*  =  {vi}f^Q  and  its  value  q*. 

Step  0.  Calculate  So{v,v',t)  for  alH  G  T  and  v,v'  eV  such  that  {v,v')  E  S 
or  do{v',t)  for  all  t  G  T  and  v'  E  V,  and  calculate  a  lower  bound  q.  Set 
i  =  0,  A(i)  =  {(no,0,  cx))}. 

Step  1.  If  )C{i)  is  empty,  replace  i  by  i  —  1.  Else,  go  to  Step  3. 

Step  2.  i  =  0,  stop:  the  last  saved  path  is  optimal  and  q  is  its  probability 
of  detection.  Else,  go  to  Step  1. 

Step  3.  Remove  from  lC{i)  the  triplet  (ut,  f,  g(nt,  f))  with  the  largest  bound 

Step  4.  If  q{vt,t)  <  Q,  go  to  Step  1.  (Current  subpath  is  fathomed.) 

Step  5.  \i  i  =  0,  replace  i  by  i  -f  1,  and  go  to  Step  3.  (/C(l)  is  populated  in 
the  network  reduction  procedure.) 

Step  6.  If  t  ^  T,  then  for  each  vertex  v  E  .F(n^),  calculate  q(vA  T  1)  from 
(II. 5)  using  bounds  do{v,  f-|-l)  or  So(vt,  +  and  add  (v,  f  +  1,  q(v,t  +  i)) 
to  /C(i  -|-  1).  Replace  i  by  i  -|-  1,  and  go  to  Step  3.  Else,  let  q  =  q{vt,  t)  and 
save  the  incumbent  path  {vi}JLq,  and  go  to  Step  1. 

We  now  present  the  network  reduction  procedure  that  can  be  imple¬ 
mented  as  part  of  Step  0  of  Algorithm  2.  The  procedure  assumes  that  a  static  bound 
do{v',  t)  is  available  as  well  as  a  lower  bound  on  the  optimal  probability  of  detection  q. 
If  the  directional  static  bound  is  available  instead  of  the  static  bound,  replace  do{v',  t) 


32 


by  So{v,v',t)  in  the  procedure  below.  We  note  again  that  the  network  reduction  pro¬ 
cedure  is  valid  under  the  assumption  that  g{vs-i,  Vs,  s)  does  not  vary  with  the  choice 
of  Us-i  £  'Ti{vs)  among  all  nodes-of-£rst-contact  {vs,s).  We  generalize  this  in  Section 

C. 


Network  Reduction  Procedure  Rl. 

Input:  A  time-expanded  network  with  arc  length  Cn,n',  the  initial  tar¬ 

get  distribution  p{-,  1),  Markov  transition  matrix  F,  and  upper  bound  q. 

Output:  Nodes-of-first-contacts  {vs,s)  which  are  not  dominated  by  others. 

Step  1.  Find  all  nodes-of-hrst-contact  {vg,  s)  G  A/”.  If  none  exist  with  s  >  1, 
then  stop. 

Step  2.  For  each  {vg,  s),  calculate  q{vs,  s)  =  p{(l>{vg),  s)g{vg^i,  Vg,  s)+do{vg,  s). 

Step  3.  Eliminate  all  nodes-of-hrst-contact  {vg,  s)  with  q{vg,  s)  <  q. 

Step  4.  For  each  nodes-of-hrst-contact  {vg,  s)  not  eliminated,  store  the  triplet 
{vg,s,q{vs,s))  in  }C{1). 

Table  2  illustrates  the  ehect  of  the  network-reduction  technique  as  ap¬ 
plied  to  the  same  problem  instances  as  in  Table  1.  Columns  5  and  6  of  Table  2 
present  the  run  time  and  number  of  branching  attempts,  respectively,  for  Algorithm 
2  with  static  bound  and  network  reduction.  On  average,  the  run  times  and  and  the 
branching  attempts  are  reduced  to  5.3%  and  5.0%  of  the  corresponding  numbers  ob¬ 
tained  without  the  network  reduction  technique  (see  columns  3  and  4  in  Table  1), 
respectively.  When  applying  both  network  reduction  and  directional  static  bound, 
we  obtain  the  run  times  and  numbers  of  branching  attempts  reported  in  columns  7 
and  8  of  Table  2.  It  is  clear  that  network  reduction  and  directional  static  bound  have 
complementary  positive  ehect  and  the  run  times  and  numbers  of  branching  attempts 
are  further  reduced. 

Table  3  presents  computational  results  for  a  larger  problem  instance 
with  15  by  15  cells  and  a  time  horizon  T  =  20.  Again,  the  searcher  starts  in  the 
upper-left  cell  and  the  targets  starts  in  the  center  cell.  As  seen  from  Table  3,  the  run 
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times  remain  rather  short  for  Algorithm  2  with  directional  static  bound  and  network 
reduction  (columns  3  and  4)  while  the  times  increases  dramatically  for  Algorithm  1 
with  dynamic  bound  (columns  5  and  6).  Furthermore,  Algorithm  2  with  directional 
static  bound  and  network  reduction  is  less  sensitive  to  the  detrimental  effect  of  a 
near  stationary  target  (e.g.,  p  =  0.9)  and  high  glimpse  detection  probability  (e.g., 
g{v,v',t)  =  0.9  or  0.99),  a  case  where  the  bounds  tend  to  be  weak. 


g{v,v',t) 

P 

Algo.  2: 

Time 

(sec.) 

D-Static  &  Red. 
Branching 
attempts 

Algo.  1: 
Time 
(sec.) 

Dynamic  Bound 
Branching 
attempts 

0.3 

3.08 

103,811 

20.24 

328,672 

0.3 

0.6 

2.94 

66,185 

21.22 

311,645 

0.9 

3.58 

238,089 

107.11 

1,352,503 

0.3 

3.14 

122,941 

20.28 

248,727 

0.6 

0.6 

2.92 

51,977 

30.47 

288,738 

0.9 

4.41 

571,649 

701.17 

8,668,034 

0.3 

3.17 

129,276 

29.61 

292,818 

0.9 

0.6 

2.95 

57,353 

45.05 

404,299 

0.9 

5.94 

1,076,948 

2323.48 

30,173,994 

0.3 

3.13 

128,776 

31.97 

301,498 

0.99 

0.6 

2.92 

60,449 

48.42 

441,664 

0.9 

5.77 

1,016,918 

3092.63 

45,329,829 

Table  3.  Run  times  and  number  of  branching  attempts  for  Algorithm  2,  with  direc¬ 
tional  static  bound  and  network  reduction,  compared  with  Algorithm  1  with  dynamic 
bound  on  15  by  15  cell  search  problem  with  time  horizon  T  =  20. 


The  same  problem  instances  were  also  examined  in  [32],  which  reports 
a  run  time  of  10.9  seconds  for  the  case  with  g{v,v',t)  =  0.6  and  p  =  0.6  using  a 
C-|--|-  implementation  of  Algorithm  1  with  a  dynamic  bound  running  on  a  2.6  GHz 
computer.  We  observe  that  our  implementation  appears  to  be  somewhat  slower  than 
the  one  achieved  in  [32]  with  30.47  seconds  compared  to  10.9  seconds  (on  a  presumedly 
slightly  slower  computer).  However,  Algorithm  2  with  directional  static  bound  and 
network  reduction  appears  to  offer  a  noticeable  advantage  over  Algorithm  1  with 
dynamic  bound  as  derived  in  [32].  In  principle.  Algorithm  1  with  dynamic  bound  can 
also  be  speeded  up  by  using  the  proposed  network  reduction  technique.  However,  we 
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have  not  examined  that  possibility.  We  adopt  a  directional  static  bound  with  network 
reduction  as  the  basis  for  extension  to  the  case  with  side  constraints  discussed  in  the 
next  section. 

C.  BRANCH-AND-BOUND  ALGORITHM  FOR  RSSP 

We  now  turn  the  attention  to  the  full  problem  with  side  constraints,  i.e.,  the 
resource-constrained  single  searcher  problem  (RSSP)  formulated  in  Section  A.  We 
first  develop  a  static  bound  based  on  Lagrangian  relaxation  that  can  be  used  within 
a  branch-and-bound  algorithm  in  the  form  of  Algorithms  1  and  2.  Second,  we  briefly 
discuss  the  development  of  a  Lagrangian  directional  static  bound.  Third,  we  develop 
a  series  of  network  reduction  techniques.  Fourth,  we  combine  the  resulting  procedures 
and  present  the  complete  algorithm. 

1.  Lagrangian  Static  Bound 

Consider  the  time-expanded  network  (A/*,  A),  see  Subsection  B.l,  with  the  arc 
rewards  Cny  given  in  (11.15).  Now,  we  also  endow  each  arc  (n,  n')  E  A,  n  =  {v,t  —  1) 
and  n'  =  with  weights  ri^n,n'  =  G  X.  While  computing  the  static 

bound  in  the  case  of  no  side  constraints  amounts  to  solving  a  longest-path  problem 
on  the  time-expanded  graph,  a  similar  bound  in  the  case  with  side  constraints  will 
need  to  account  for  those  constraints.  Specihcally,  a  static  bound  can  be  obtained  by 
solving  a  constrained  longest-path  problem.  We  formulate  this  problem  as  an  integer 
program  on  a  slightly  modihed  time-expanded  graph. 

We  use  the  same  time-expanded  graph  as  in  Subsection  B.l,  with  the  following 
modihcations.  We  recall  that  V  is  the  set  of  vertices  where  the  search  can  end.  Now, 
every  node  n  =  (v,  t),  w  G  V,  t  G  T  is  connected  to  the  artificial  destination  node  h 
with  an  arc.  This  modihcation  allows  the  searcher  to  terminate  the  search  prior  to 
time  period  T  to  avoid  violating  the  side  constraints.  Furthermore,  all  arcs  (n,  n)  with 
n  =  {v,T),  V  ^  V,  are  removed  from  the  time-expanded  network.  This  modihcation 
makes  the  searcher  end  at  w  G  V.  We  still  let  A  denote  the  set  of  all  arcs. 
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We  formulate  the  constrained  longest-path  problem  on  the  time-expanded 
graph  (A/*,  A)  as  an  integer  program.  We  consider  an  ordering  of  A  and  let  A  denote 
the  lA/”!  by  |Al|  node-arc  incidence  matrix  for  the  time-expanded  graph.  For  each  arc 
a  =  {n,n')  G  A,  let  =  1,  A„',a  =  —1,  and  =  0  for  any  n"  G  Af,n"  ^  n,  n' 

be  the  elements  of  A.  Let  b  denote  the  |A/’|-vector  with  =  lAh  =  —1  and  6^  =  0 
for  all  n  G  A/'\{no,h}.  We  also  dehne  the  additional  notation;  f  =  (ri,r2, . . .  i 
where  T  denotes  the  transpose.  We  collect  the  rewards  Cn,n'  ia  the  | A | -dimensional 
row  vector  c.  Also,  for  each  i  G  X,  we  define  the  |A|-dimensional  row  vector  r*  to 
contain  the  weights  ri^n,n'  and  we  let  R  be  the  |X|  by  |A|  matrix  with  as  its  rows. 
Finally,  we  let  x  be  a  |A|-dimensional  column  vector,  where  Xn,n'  =  1  if  arc  {n,n')  is 
used  by  a  path,  and  zero  otherwise.  Then,  the  constrained  longest-path  problem,  see 
[1],  can  be  written  as: 


=  max  cx 

(11,18) 

S.t.  Ax  =  b 

Rx  <  r. 

(11.19) 

In  principle,  the  solution  of  the  constrained  longest-path  problem  provides  a  static 
bound.  However,  the  constrained  longest-path  problem  is  NP-complete  even  for  the 
case  with  an  acyclic  graph  ([22],  pages  213-214).  Hence,  we  prefer  to  avoid  solving  such 
problems  within  an  algorithm.  We  proceed  by  introducing  an  additional  relaxation 
that  is  motivated  by  the  solution  approach  for  the  constrained  shortest-path  problem 
in  [24,  9,  10]. 

Using  the  standard  theory  of  Lagrangian  relaxation  (see,  e.g.,  [1],  Chapter  16) 
and  an  |X|-dimensional  row  vector  A  >  0,  we  hnd  that 

<  ^(A)  =  max  cx  —  A(Rx  —  f)  (11.20) 

xe{o,i}l-^ 

s.t.  Ax  =  b. 
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Rewriting  the  objective  function,  we  can  optimize  the  Lagrangian  upper  bound  2;(A) 
through 


2;  =  min2:(A)  (11.21) 

A>0 

=  min  max  (c  — AR)x  +  Af  (11.22) 

A>o  xe{o,i}i-^ 

s.t.  Ax  =  b. 

For  any  hxed  A  >  0,  computing  the  upper  bound  ;2(A)  simply  requires  the  solution  of 
a  longest-path  problem  with  Lagrangian-modihed  arc  lengths  in  an  acyclic  graph.  The 
outer  minimization  over  A  can  be  solved  by  several  methods  [19,  5,  14,  9,  10].  Since 
we  anticipate  only  a  small  number  of  side  constraints,  it  suffices  to  use  a  repeated 
coordinate  search.  Given  an  optimal  or  near-optimal  A,  we  obtain  the  Lagrangian 
static  bound  by  carrying  out  one  backward  longest-path  calculation  in  the  time- 
expanded  graph  from  node  n  to  all  nodes  n  &  M  using  the  Lagrangian  modified 
arc  reward  c  —  AR.  More  specihcally,  an  arc  {n,n')  =  {{v,s),{v',s  +  1))  G  A, 
s  =  1,  2, ...,  T  —  1,  is  endowed  with  the  reward 

Cn,n'  ^  i,n,n' ;  (11.23) 

i£l 

where  Cn,n'  is  given  by  (11.15). 

We  also  derive  and  implement  a  Lagrangian  directional  static  bound  similar  to 
the  one  in  Subsection  B.2c.  However,  the  derivation  is  a  straightforward  combination 
of  Subsection  B.2c  and  the  approach  described  above.  Hence,  we  omit  it.  This 
derivation  results  in  a  similar  Lagrangian  problem  to  the  one  in  (11.21),  but  now 
defined  on  the  node-and-time  expanded  network.  We  still  refer  to  the  Lagrangian 
multiplier  as  A  and  the  Lagrangian  upper  bound  as  2: (A).  We  denote  the  Lagrangian 
directional  static  bound  computed  in  this  way  by  6o{v,v',t)-  Since  the  Lagrangian 
directional  static  bound  is  at  least  as  strong  as  the  Lagrangian  static  bound,  we  carry 
out  the  Lagrangian  relaxation  only  in  the  node-and-time  expanded  network  to  hnd  a 
Lagrangian  multiplier  A  >  0. 
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2.  Network  Reduction 

We  propose  and  examine  three  techniques  for  reducing  the  size  of  the  network 
prior  to  application  of  a  branch-and-bound  procedure.  First,  we  use  dominance  rules 
to  eliminate  edges  that  cannot  be  on  an  optimal  path.  Second,  we  describe  the 
application  of  “preprocessing”  techniques  frequently  used  prior  to  solving  constrained 
shortest-path  problems.  Third,  we  modify  the  procedure  described  in  Subsection 
B.2d. 

a.  Edge  Dominance 

There  are  several  situation  where  a  vertex  v'  G  Elv)  can  be  eliminated 
as  candidate  for  visit  from  vertex  v.  Such  “dominance  tests”  are  case  dependent,  but 
can  be  effective  in  reducing  the  number  of  edges.  We  describe  one  situation  where  we 
use  “edge  dominance.” 

In  many  practical  situations,  there  are  two  resources:  risk  and  fuel. 
If  higher  altitude  implies  lower  risk  and  lower  glimpse  detection  probability,  and 
climbing  to  higher  altitude  consumes  more  fuel  than  level  flight,  then  we  can  eliminate 
some  edges  in  the  graph  (V, T)  by  edge  dominance.  Suppose  that  fi{v,v',t)  and 
/2(n,  v',  t)  represent  risk  and  fuel,  respectively.  Also  suppose  that  the  risk  /i(n,  n',  t)  = 
/i(n'),  i.e.,  only  depends  on  v' .  Let  'ipiv)  be  the  altitude  of  waypoint  n  G  V.  Then,  if 
we  have  the  above  described  situation,  we  use  the  following  (one-step)  procedure  to 
reduce  the  size  of  the  graph  (V,T). 

Edge  Dominance  Procedure  R2. 

Step  1.  Delete  any  edge  (n,  v')  G  £  that  satisfies  fi{v')  =  0  and  'ipiy)  <  ipiy'). 

We  note  that  Procedure  R2  takes  advantage  of  the  fact  that  if  there  is  no  risk  at  n', 
then  there  is  no  need  to  increase  altitude  when  moving  from  v  to  v' .  The  altitude 
can  be  increased  later  if  need  be. 
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b.  Preprocessing 

It  is  well  known  that  the  side  constraints  (11.19)  can  be  used,  prior 
to  any  main  calculations,  to  identify  nodes  and  arcs  in  the  time-expanded  network 
{Af,A)  that  cannot  lie  on  any  feasible  path  [4,  15,  9,  10].  Such  nodes  and  arcs 
can  be  eliminated  from  which  reduces  the  size  of  the  problem  that  needs  to 

be  considered  when  computing  bounds  and  carrying  out  branching.  Moreover,  the 
reduction  of  the  time-expanded  network  typically  also  strengthens  the  Lagrangian 
relaxation,  i.e.,  reduces  the  gap  z  —  z*,  (see  (11.21)  and  (11.20)),  and,  hence,  reduces 
the  need  for  branching.  We  adopt  the  follow  procedure,  adapted  from  [9,  10],  to  carry 
out  arc  preprocessing: 

Preprocessing  Procedure  R3. 

Input:  A  time-expanded  network  {Af,A)  and  arc  lengths 

Output:  The  updated  (or  reduced)  time-expanded  network  {J\f,A). 

Step  1.  Set  number  of  iterations  k  and  k  =  1. 

Step  2.  For  alH  G  X  and  n  G  A/”,  compute  a  minimum-weight  rio-n  subpath 
distance  Riin)  and  a  minimum- weight  n-fi  subpath  distance  rj(n)  in  (A/*,  A) 
with  respect  to  weights 

Step  3.  Delete  any  arc  (n,  n')  G  A  with 

Ri{n)  +  +  ri{n')  >  R  for  any  i  G  X.  (11.24) 

Step  4.  If  k  <  k  and  at  least  one  arc  was  deleted  in  Step  3,  replace  khy  k  +  1, 
and  go  to  Step  2.  Else,  stop. 

If  a  lower  bound  on  the  probability  of  detection  is  available,  we  also 
carry  out  similar  preprocessing  with  respect  to  arc  reward  Cny  and  the  Lagrangian 
modihed  arc  reward  Cny  =  Cny  —  ^i^i,ny ,  see  [9,  10]. 

We  describe  the  preprocessing  procedure  for  the  time-expanded  net¬ 
work  and  argue  that  it  improves  the  Lagrangian  static  bound.  However,  the  same 
methodology  applies  to  the  node-and-time  expanded  network  and  it  improves  the 
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Lagrangian  directional  static  bound.  In  our  main  algorithm  (described  in  Subsection 
C.3),  we  also  apply  preprocessing  to  the  node-and-time  expanded  network  and  denote 
that  procedure  R3'. 


c.  Vertex  Dominance  for  Distant  Target 

As  in  Subsection  B.2d,  we  consider  the  case  where  the  searcher’s  and 
the  target’s  locations  are  initially  some  distance  apart  and  derive  a  network  reduction 
procedure  that  takes  advantage  of  that  situation.  In  contrast  to  Subsection  B.2d, 
the  subpath  used  to  reach  a  node-of-£rst-contact  is  now  important  since  the  resource 
consumption  along  different  subpaths  may  be  different.  Hence,  (uq,  0)  now  needs  to  be 
connected  to  each  node-of-hrst-contact  {vg,  s)  with  multiple  “jump”  arcs  representing 
the  different  possible  subpaths  and  resource  consumptions  used  to  reach  {vs,s).  A 
standard  path  enumeration  algorithm  (see,  e.g.,  [11])  can  enumerate  the  different 
subpaths,  at  least  as  long  as  s  is  relatively  small.  Multiple  arcs  to  a  node-of-hrst- 
contact  can  also  be  used  to  model  the  situation  with  edge-dependent  glimpse  detection 
probability,  a  case  ignored  in  Subsection  B.2d. 

After  all  the  arcs  are  generated  to  all  the  nodes-of-£rst-contact,  a  num¬ 
ber  of  them  can  be  deemed  uninteresting  and  can  be  eliminated  using  dominance  rules 
of  the  form:  If  an  arc  from  {vq,  0)  to  {Vg,  s)  has  no  larger  probability  of  detection  and 
no  smaller  consumption  of  each  resource  as  another  parallel  arc  and  the  two  arcs  are 
not  identical,  the  hrst  arc  is  dominated  and  can  be  eliminated.  In  sets  of  identical 
parallel  arcs,  we  also  eliminate  all  but  one.  Trivially,  arcs  with  resource  consumption 
greater  than  the  specified  limits  are  also  removed.  Moreover,  if  a  lower  bound  on 
the  optimal  probability  of  detection  exists,  it  can  be  used  to  eliminate  more  arcs  as 
described  in  the  following. 

Below  we  formally  describe  this  network  reduction  procedure  based  on 
vertex  dominance  for  a  distant  target.  We  note  that  the  procedure  is  more  effective 
after  (i)  applying  network  reduction  procedures  R2,  R3  and  R3',  (ii)  hnding  A  that 
(approximately)  optimizes  the  Lagrangian  upper  bound  ;2(A)  on  the  node-and-time 
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expanded  network,  and  (iii)  computing  the  directional  static  bound  So{v,v',t)  and 
the  Lagrangian  directional  static  bound  So{v,v',t)  for  each  node  {v',t)  G  J\f.  So 
we  assume  that  these  calculations  have  been  carried  out.  During  these  calculations, 
feasible  paths  may  be  obtained.  Such  paths  provide  lower  bounds  on  the  optimal 
probability  of  detection.  Let  q  denote  the  largest  lower  bound  found  so  far. 

Vertex  Dominance  Procedure  R4. 

Input:  A  time-expanded  network  arc  lengths  c„^„/  and  Tn^n'- 

Output:  The  nodes-of-£rst-contact  {vg,  s)  which  are  not  dominated  by  others. 

Step  1.  Find  all  nodes-of-hrst-contact  {vs,s)  G  A/”.  If  none  exist  with  s  >  1, 
then  stop. 

Step  2.  For  each  node-of-first-contact  {vs,s),  enumerate  all  subpaths  (no,0) 
to  {vs,s). 

Step  3.  For  each  node-of-first-contact  {vs,s)  and  subpath  V  =  carry 

out  the  tests: 

If  any  of  the  following  is  true,  then  eliminate  V: 

(i)  For  some  remaining  {vo,0)-{vs,  s)  subpath  V',  ViiV)  >  rj('P')  for  all 
i  G  X  and  q{V)  <  q{V'). 

(ii)  q{V)  +  5o{vs-i,  Vg,  s)  <  q. 

(iii)  r-iiV)  +  s))  >  fi  for  some  i  G  X. 

Step  3  can  also  be  augmented  with  a  test  on  the  Lagrangian-modified  probability  of 
detection  using  6o{v,v',t)  if  a  near-optimal  multiplier  A  is  available. 

3.  Algorithm 

We  now  state  the  complete  algorithm  for  RSSP.  The  algorithm  starts  with 
network  reductions  procedures  R2,  R3,  and  R3'.  The  next  step  solves  the  Lagrangian 
problem  (11.22)  and  determines  a  near-optimal  A.  (The  calculations  are  actually 
carried  out  in  the  node-and-time  expanded  network  as  we  prefer  to  use  the  Lagrangian 
directional  static  bound.)  If  a  feasible  path  becomes  available  during  the  procedures 
described  above,  network  reduction  procedure  R3'  is  repeated  now  using  checks  with 
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respect  to  arc  reward  and  its  Lagrangian  modified  arc  reward.  The  next  steps 
are  to  compute  the  directional  static  bound  of  Subsection  B.2c  (i.e.,  So{v,v' ,t)),  the 
Lagrangian  directional  static  bound  as  described  in  Subsection  C.l  (i.e.,  6o{v,v' ,t)), 
and  bounds  on  resource  consumption  along  any  path  extension  (i.e.,  rj(n),  i  G  X).  The 
final  steps  before  the  branch-and-bound  procedure  is  to  implement  network  reduction 
procedure  R4. 

We  implement  the  branch-and-bound  procedure  as  an  implicit  path-enumeration 
in  the  time-expanded  network.  The  procedure  amounts  to  a  depth-first  search  cou¬ 
pled  with  optimality  and  feasibility  checks  using  the  computed  bounds,  which  follows 
[9,  10].  The  complete  algorithm  takes  the  following  form: 

Algorithm  3  (Solves  RSSP). 

Input:  A  time-expanded  network  (A/*,  A)  with  arc  lengths  Cn^n'  and  Vny  and  a 
node-and-time  expanded  network  (S,  hi)  with  arc  lengths  and  The 
initial  target  distribution  p(-,l),  Markov  transition  matrix  T,  and  upper 
bound  q. 

Output:  An  optimal  search  path  V*  =  {vi}f^Q  and  its  value  q*. 

Step  1.  Apply  network  reduction  procedures  R2,  R3  and  R3'. 

Step  2.  Find  A  that  approximately  optimizes  the  Lagrangian  upper  bound 
;2(A)  in  the  node-and-time  expanded  graph  (S,  hi).  If  a  feasible  solution  is 
found,  set  the  probability  of  detection  on  the  corresponding  path  equal  to 
q.  Otherwise,  set  q  =  —oo. 

Step  3.  If  a  feasible  solution  is  found  so  far,  implement  R3'  also  with  respect 
to  arc  reward  and  Lagrangian  modified  arc  reward  using  q. 

Step  4.  Ignoring  side  constraints,  compute  the  directional  static  bound  50(1^,  v',  t) 
for  all  expanded  nodes  ^  =  (n,  n')  in  (S,  O),  with  n  =  (n,  t  —  1),  n'  =  (n',  t), 
n,  v'  e  V. 

Step  5.  Using  A  from  Step  2,  compute  the  Lagrangian  directional  static 
bound  (5o(n,U,t)  for  all  expanded  nodes  ^  =  {n,n')  in  (S,r2),  with  n  = 

{v,  t  —  1),  n'  =  {v',  t),  V,  v'  e  V. 

Step  6.  For  each  i  G  X,  compute  the  minimum  distance  rj(n)  from  each  node 
n  &  M  back  to  n  by  solving  a  single,  backwards,  shortest-path  problem  in 
the  time-expanded  graph  {J\f,A)  starting  from  h  using  arc  lengths  ri^n,n'- 
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Step  7.  Apply  network  reduction  procedure  R4. 

Step  8.  Apply  a  standard  path-enumeration  procedure  (e.g.,  [11])  in  {Af^A) 
with  the  following  modihcations: 

(i)  The  path-enumeration  commences  from  no,  but  extends  a  current  sub¬ 
path  {niYi^Q  along  arc  {nt,n)  =  ((nt,t),(n,t  -|-  1))  if  and  only  if  the 
following  conditions  hold; 

•  For  all  i  E  I,  {no,  ni, ...,  nt,  n}  can  be  extended  to  a  path  whose  i-th 
resource  does  not  exceed  ri,  i.e., 

ri({no,  ni, ...,  rit,  n})  rj(n)  <  A.  (11.25) 

•  (no,  rii,  can  be  extended  to  a  path  with  probability  of  de¬ 

tection  exceeding  q,  i.e.. 


g({no,  ni, ...,  n*,  n})  -h  v,t  +  l)  >  q-  (11.26) 

•  (no,  rii,  can  be  extended  to  a  path  whose  Lagrangian- modified 

probability  is  no  less  than  q,  i.e., 

t 

g({no,ni,  ...,nt,n})  - 

iex  1=1 

-  E  +  Ar  5o(^'i,  W,  t  +  1)  >  ^  (11.27) 

i£l 

(ii)  Whenever  the  algorithm  identifies  a  path  V  with  q{V)  >  q  and  ri{V)  < 
fi,i  G  X,  replace  q  by  q(V). 


In  Step  8,  the  checks  (11.25),  (11.26),  and  (11.27)  prevent  the  enumeration 
of  paths  that  can  be  determined,  using  the  computed  bounds,  to  not  be  optimal. 
Specifically,  (11.25)  prevents  the  extension  of  subpaths  that  cannot  result  in  a  feasible 
path  with  respect  to  the  side  constraints.  Since  v,t  +  Y  is  a  valid  upper  bound 
on  the  probability  of  detection  during  time  periods  t  -|-  2,  t  -|-  3,  ...,  T,  the  left-hand 
side  of  (11.26)  is  an  upper  bound  on  the  probability  of  detection  along  any  path  that 
starts  with  the  subpath  (no,  ni, ...,  n^,  n}.  Hence,  the  subpath  cannot  be  extended 
to  a  path  with  larger  probability  of  detection  than  q  if  (11.26)  fails.  In  (11.27),  the 
probability  of  detection  along  {no,ni,  is  modihed  by  Lagrangian  terms  and 
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the  Lagrangian  directional  static  bound  is  used.  The  resulting  check  can  be  shown 
to  be  valid  using  standard  Lagrangian  relaxation  theory  and  the  argument  above. 

We  also  implement  Step  8  with  a  “branching  strategy”  based  on  the  La¬ 
grangian  directional  static  bound.  Specihcally,  we  first  consider  extending  the  current 
subpath  {riiYi^Q  along  the  arc  {nt,n)  with  the  largest  Lagrangian  directional  static 
bound  among  all  the  nodes  in  the  forward  star  of  rit.  Second,  we  consider  extend¬ 
ing  {riiW^Q  with  the  node  corresponding  to  the  second  largest  Lagrangian  directional 
static  bound,  etc.  This  branching  strategy  is  analogous  to  the  one  in  Step  3  of  Algo¬ 
rithms  1  and  2.  We  have  also  experimented  with  using  the  directional  static  bound 
instead  of  the  Lagrangian  directional  static  bound  and  found  it  to  usually  result  in 
comparable  run  times.  However,  the  Lagrangian  directional  static  bound  appears 
faster,  on  average. 

Since  Algorithm  3,  in  the  worst  case,  enumerates  all  feasible  paths,  it  is  guar¬ 
anteed  to  hnd  an  optimal  solution  of  RSSP. 

D.  NUMERICAL  EXAMPLE 

This  section  describes  computational  experiments  with  Algorithm  3  applied 
to  RSSP  with  side  constraints  on  risk  exposure  and  fuel  consumption.  We  carry  out 
all  experiments  on  the  same  computational  platform  as  in  Section  B. 

We  consider  a  military  planning  situation  where  a  UAV  is  assigned  a  mission 
to  search  and  detect  a  high-value  moving  target.  Planners  wish  to  determine  a  flight 
path  over  the  area  of  interest  (AOI)  that  maximizes  the  probability  of  detecting  the 
target.  The  UAV  will  start  its  path  at  a  known  waypoint  with  a  known  fuel  supply, 
and  will  return  to  the  same  waypoint  before  the  fuel  tank  is  empty.  Doctrine  specihes 
that  the  UAV  cannot  be  assigned  a  path  with  higher  risk  than  a  specihc  threshold. 
The  AOI  is  partially  under  enemy  control  and  any  aircraft  flying  over  the  AOI  could  be 
shot  down  by  enemy  surface-to-air  missiles  (SAMs),  anti-aircraft  artillery,  and  small- 
arms  hre.  Flying  at  a  high  altitude  would  reduces  that  risk,  but  it  will  also  reduce  the 
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quality  of  UAV’s  sensor  readings.  Consequently,  the  UAV  may  change  altitude  during 
the  course  of  the  mission  to  balance  risk  and  sensor  quality.  Changing  from  low  to 
high  altitude  consumes  more  fuel  than  level  flight.  Hence,  the  number  of  time  periods 
available  for  search  depends  on  the  fuel  consumption  and  therefore  the  vertical  flight 
prohle. 

We  model  this  situation  by  dividing  the  AOI  into  10  by  10  cells  (see  Figure 
6).  The  airspace  over  each  cell  is  vertically  discretized  into  two  altitudes  (’’low  and 
high”).  The  heavily  shaded  cells  (Ci)  in  Figure  6  represent  an  urban  area  over  which 
the  UAV’s  risk  is  high  and  its  glimpse  detection  probability  is  low.  The  unshaded  cells 
(C3)  represent  open  terrain  where  there  is  no  risk  and  the  UAV’s  glimpse  detection 
probability  is  high.  The  lightly  shaded  cells  (C2)  represent  an  area  with  intermediate 
risk  and  intermediate  glimpse  detection  probability.  We  note  that,  while  this  problem 
instance  does  not  directly  relate  to  any  real-world  operational  situation,  we  believe 
that  it  illustrate  a  fairly  practical  situation. 

We  assume  that  the  risks  at  different  edges  along  a  path  are  independent.  If 
the  probability  of  the  UAV  surviving  edge  (n,  v')  G  ^  is  a(v,  v'),  then  the  probability 
of  surviving  the  path  is  simply  nfLia(vi-i,  vi).  Let  d  be  a  lower  limit  on 

the  survival  probability.  Then,  a  standard  logarithmic  transformation  leads  to  the 
following  constraint 

k 

-  log  Vi)  <  -  log  d,  (11.28) 

i=i 

which  is  in  the  form  (II. 3)  and  (II.4)  with  fi{vi-i,vi,l)  =  —\oga{vi-i,vi)  and  A  = 
-  log  d. 

For  this  compntational  experiment,  we  assume  that  the  glimpse  detection 
probability  and  the  snrvival  probability  for  an  edge  (n,  v')  G  S  depend  only  on  the 
cell  and  altitnde  corresponding  to  vertex  n'  G  V  as  listed  in  Table  4.  We  note  that 
glimpse  detection  probability  at  high  altitude  is  assumed  to  be  70%  of  that  at  low 
altitude  and  the  failnre  probability  (complement  of  the  survival  probability)  at  high 
altitude  is  30%  of  that  at  low  altitude. 
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Figure  6.  An  area  of  interest  composed  of  10  by  10  cells  and  two  altitudes.  Heavily 
shaded  cells  (Ci),  lightly  shaded  cells  (C2),  and  unshaded  cells  (C3)  describe  risky, 
moderately  risky,  and  non-risky  area,  respectively.  The  circle  indicates  the  cell  over 
which  searcher  starts,  and  the  triangle  specihes  the  initial  position  of  the  target. 

The  UAV  enters  the  airspace  at  high  altitude  over  the  northwest  cell  (cell  1; 
cells  are  numbered  from  left  to  right,  and  from  top  to  bottom)  and  will  return  to 
the  same  cell  at  either  high  or  low  altitude  at  the  end  of  the  mission.  The  searcher 
is  located  at  one  vertex  each  time  period  and  searches  the  corresponding  cell.  For 
the  next  time  period,  the  searcher  can  stay  at  the  same  vertex,  change  altitude  over 
the  same  cell,  or  move  to  a  vertex  (at  any  altitude)  corresponding  to  a  vertically  or 
horizontally  adjacent  cell.  The  maximum  number  of  time  periods  is  T  =  40,  but  the 
fuel-consumption  constraint  may  limit  the  number  of  periods  to  less  than  40.  We 
assume  that  the  fuel  consumption  at  each  time  step  is  as  follows:  10  units  if  there 
is  no  altitude  change,  12(9)  units  if  changing  from  low(high)  altitude  to  high(low) 
altitude.  The  initial  position  of  the  target  is  the  center  of  the  high  risk  region  (cell 
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Cells 

Altitude 

Glimpse  probability 

Survival  probability 

Cl 

low 

0.20 

0.960 

high 

0.14 

0.988 

C2 

low 

0.40 

0.980 

high 

0.28 

0.994 

Cs 

low 

0.60 

1.000 

high 

0.42 

1.000 

Table  4.  Glimpse  detection  probability  g{v,v',t)  and  survival  probability  a{v,v')  for 
different  cells  and  altitude. 


68).  The  target  remains  in  the  current  cell  with  a  probability  p  =  0.6  for  the  next 
time  period  or  moves  to  one  of  the  vertically  or  horizontally  adjacent  cells  with  equal 
probability. 

The  survival-probability  limit  is  a  threshold  that  is  set  by  the  commander  or 
planner.  A  search  path  with  lower  survival  probability  than  the  threshold  would  not 
be  accepted.  In  this  experiment,  we  consider  the  survival  probability  limits  0.95, 
0.90,.. .,0.70,  and  fuel  consumption  limit  300,  325,  ...,  450. 

We  solve  this  problem  instance  using  Algorithm  3.  Tables  5,  6  and  7  report 
computational  results  for  different  combinations  of  survival  probability  and  fuel  limits 
for  the  UAV.  When  the  fuel  limit  is  tight  (e.g.,  300  and  325),  the  UAV  cannot  operate 
for  the  full  duration  of  40  time  steps.  We  observe  that  increasing  the  fuel  limit  beyond 
425  does  not  increase  the  probability  of  detection  as  the  time  limit  of  40  periods 
becomes  active.  The  average  run  time  is  580  seconds,  with  a  standard  deviation  of 
792.  All  problem  instances  are  solved  within  one  hour  and  typically  in  much  less. 
Figure  7  shows  the  optimal  path  given  survival  probability  limit  0.90  and  fuel  limit 
400.  The  solid  lines  and  the  dashed  lines  represent  flight  segment  at  low  and  high 
altitude,  respectively. 

We  also  consider  the  case  with  edge-dependent  glimpse  detection  probability. 
Consider  the  same  situation  as  earlier  described,  but  now  assume  that  a  move  to  a  new 
waypoint  results  in  a  lower  glimpse  detection  probability  than  if  the  searcher  already 
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was  at  that  waypoint.  Specifically,  if  n  =  v',  we  let  the  glimpse  detection  probability 
g{v,v',t)  be  as  in  Table  4;  otherwise  we  replace  g{v,v',t)  by  0.1g{v,v' ,t).  Figure  8 
shows  an  optimal  path  found  given  survival  probability  and  fuel  limits  of  0.90  and 
400,  respectively.  In  contrast  to  the  case  with  edge-independent  glimpse  detection 
probability  (Figure  7),  the  searcher  now  tends  to  stay  for  multiple  time  periods  at  the 
same  waypoints  in  high-probability  regions  to  reap  the  benehts  of  the  corresponding 
high  glimpse  detection  probability.  The  run  times  (not  reported  in  detail)  for  the 
case  with  edge-dependent  glimpse  detection  probabilities  are,  on  average,  53  seconds, 
with  a  standard  deviation  of  71  seconds.  The  reduction  in  run  time  compared  to 
the  edge-independent  case  is  caused  by  the  often  lower  glimpse  detection  probability 
(0.l5f(n,  F,  t)),  which  tightens  the  bound. 


Fuel 

limit 

Survival  prob.  limit  = 
Prob.  Survival  Fuel 

Detection  Prob. 

0.95 

Run  time 
(sec.) 

Fuel 

limit 

Survival  prob.  limit  = 
Prob.  Survival  Fuel 

Detection  Prob. 

0.90 

Run  time 
(sec.) 

300 

0.131501 

0.962466 

300 

26.89 

300 

0.135098 

0.915157 

300 

29.58 

325 

0.153228 

0.952892 

321 

98.20 

325 

0.166933 

0.901925 

322 

32.66 

350 

0.176511 

0.952892 

350 

3313.64 

350 

0.194440 

0.915157 

350 

115.94 

375 

0.204686 

0.952892 

371 

661.43 

375 

0.217277 

0.901925 

372 

84.00 

400 

0.236651 

0.962466 

400 

664.96 

400 

0.246767 

0.902268 

400 

537.71 

425 

0.237081 

0.952892 

401 

2721.68 

425 

0.249996 

0.901925 

402 

361.60 

450 

0.237081 

0.952892 

401 

2724.22 

450 

0.249996 

0.901925 

402 

361.52 

Table  5.  Computational  results  (probability  of  detection,  survival  probability,  fuel 
consumption  and  run  times)  for  Algorithm  3  with  survival  probability  limit  =  0.95 
and  0.90. 
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Survival  prob.  limit  = 
Prob.  Survival  Fuel 

Detection  Prob. 

0.85 

Run  time 
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Fuel 

limit 

Survival  prob.  limit  = 
Prob.  Survival  Fuel 

Detection  Prob. 

0.80 

Run  time 
(sec.) 

300 

0.145309 

0.851528 

300 

25.31 

300 

0.145809 

0.805953 

299 

22.03 

325 

0.173000 

0.851528 

320 

37.81 

325 

0.173218 

0.839535 

319 

37.39 

350 

0.203183 

0.851528 

350 

67.69 

350 

0.204891 

0.805953 

349 

64.17 

375 

0.222393 

0.851528 

370 

115.24 

375 

0.223377 

0.805953 

369 

116.69 

400 

0.255481 

0.851528 

400 

510.68 

400 

0.255813 

0.827710 

399 

880.65 

425 

0.255481 

0.851528 

400 

508.50 

425 

0.255813 

0.827710 

399 

880.40 

450 

0.255481 

0.851528 

400 

508.96 

450 

0.255813 

0.827710 

399 

880.46 

Table  6.  Computational  results  (probability  of  detection,  survival  probability,  fuel 
consumption  and  run  times)  for  Algorithm  3  with  survival  probability  limit  =  0.85 
and  0.80. 


Fuel 

limit 

Survival  prob.  limit  = 
Prob.  Survival  Fuel 

Detection  Prob. 

0.75 

Run  time 
(sec.) 

Fuel 

limit 

Survival  prob.  limit  = 
Prob.  Survival  Fuel 

Detection  Prob. 

0.70 

Run  time 
(sec.) 

300 

0.150092 

0.753377 

300 

21.33 

300 

0.150277 

0.742767 

299 

21.38 

325 

0.175850 

0.753377 

320 

37.50 

325 

0.176068 

0.742767 

319 

31.16 

350 

0.205935 

0.753377 

350 

63.84 

350 

0.206200 

0.742767 

349 

53.42 

375 

0.223703 

0.753377 

370 

125.97 

375 

0.224003 

0.713056 

369 

121.33 

400 

0.255813 

0.827710 

399 

1220.30 

400 

0.255813 

0.827710 

399 

1280.52 

425 

0.255813 

0.827710 

399 

1217.21 

425 

0.255813 

0.827710 

399 

1277.52 

450 

0.255813 

0.827710 

399 

1218.91 

450 

0.255813 

0.827710 

399 

1282.42 

Table  7.  Computational  results 
consumption  and  run  times)  for 
and  0.70. 


(probability  of  detection,  survival  probability,  fuel 
Algorithm  3  with  survival  probability  limit  =  0.75 
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Figure  7.  An  optimal  path  for  survival  probability  limit  0.90  and  fuel  limit  400.  The 
solid  lines  and  the  dashed  lines  represent  flight  segments  at  low  and  high  altitude, 
respectively.  (See  Figure  6  for  a  description  of  the  underlying  figure.) 
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Figure  8.  An  optimal  path  for  a  case  with  edge-dependent  glimpse  probability,  sur¬ 
vival  probability  limit  0.90,  and  fuel  limit  400.  The  solid  lines  and  the  dashed  lines 
represent  flight  segments  at  low  and  high  altitude,  respectively.  (See  Figure  6  for  a 
description  of  the  underlying  figure.) 
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III.  MULTIPLE-SEARCHER  PROBLEM 

This  chapter  addresses  the  optimization  of  multiple  searchers  with  path-  and 
time-constraints  (MSP).  We  aim  to  determine  a  path  for  each  searcher  that  maximize 
the  probability  of  detecting  a  moving  target  within  a  mission  time  by  at  least  one 
searcher.  We  refer  to  such  a  set  of  paths  as  a  “search  plan.”  We  start  by  describing 
and  formulating  MSP.  Then,  we  present  three  algorithms  (an  exact  procedure  and 
two  heuristics)  to  hnd  an  optimal,  near-optimal,  or  “good”  search  plan.  The  chapter 
also  includes  computational  studies. 

A.  PROBLEM  DESCRIPTION  AND  FORMULATION 

The  multiple-searcher  problem  (MSP)  is  an  extension  of  SP  which  considers 
a  hnite  set  of  searchers  J  =  {1,2,  ...,J}.  In  this  section,  for  ease  of  reference,  we 
formulate  MSP  by  retracing  the  relevant  portions  of  Section  A  in  Chapter  II. 

The  area  of  interest  is  discretized  into  a  hnite  set  of  cells  C  =  (1, . . . ,  C*}  and  the 
time  horizon  is  discretized  into  a  hnite  set  of  time  periods  T  =  {1,2,  ...,T}.  A  target 
occupies  one  cell  each  time  period  and  moves  among  cells  according  to  a  Markov 
process  with  known  transition  matrix  P.  Let  p{-,t)  =  \p{l,t),p{2,t), . . .  ,p{C,t)], 
where  p{c,  t)  is  the  probability  that  the  target  is  in  cell  c  G  C  at  the  beginning  of  time 
period  t  &  T  and  the  target  has  not  been  detected  before  t  by  any  searcher.  We  refer 
to  p(-,t)  as  the  undetected  target  distribution.  The  initial  target  distribution  p(-,  1) 
is  known. 

The  multiple  searchers  move  through  a  designated  airspace  over  the  area  of 
interest  with  the  goal  of  hnding  the  moving  target  on  the  ground.  The  airspace  over 
each  cell  c  G  C  is  vertically  discretized  into  a  set  of  altitudes  Ti  =  {1,2,...,  H}.  For 
any  c  G  C  and  h  eH,  we  refer  to  the  cell-altitude  pair  (c,  h)  as  a  waypoint  where  the 
searchers  can  loiter  and  carry  out  search  of  cell  c.  We  model  the  designated  airspace 
by  a  directed  network  (V,T),  with  set  of  vertices  V  and  set  of  directed  edges  S,  in 
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which  vertices  v  =  {c,  h)  E  V  represent  waypoints  and  directed  edges  e  =  {v,  v')  G  £ 
represent  transition  between  waypoints  n,  v'  G  V.  The  searchers  can  only  transit 
between  two  waypoints  that  are  located  physical  adjacent  to  each  other.  Let  J^{v)  dV 
be  the  set  of  vertices  that  are  adjacent  to  n  G  V.  We  adopt  the  convention  that 
V  G  J^{v)  for  all  n  G  V.  Then,  the  set  of  edges  £  =  {(v,n')  G  V  x  V|  v'  G 

During  each  time  period  t  E  T,  each  searcher  is  located  at  a  particular  vertex 
(waypoint).  Any  number  of  searchers  can  occupy  a  vertex  in  the  same  time  period. 
We  also  assume  that  there  is  no  transit  time  between  waypoints.  Hence,  (n,  v')  E  £ 
simply  represents  search  at  waypoint  v  followed  by  search  at  waypoint  v'  in  the  next 
time  period.  We  note  that  the  edge  (n,  v)  E  £  represents  searching  at  waypoint  v  for 
two  consecutive  time  periods. 

Let  0  ;  V  — >  C  be  the  function  that  specifies  the  cell  over  which  a  vertex 
is  located,  i.e.,  cell  0(v)  is  searched  from  vertex  v.  We  denote  the  initial  vertex 
(waypoint)  of  searcher  j  E  J  prior  to  time  period  1  by  Ug  ^  V.  We  also  define  W  C  V 
to  be  a  set  of  possible  destination  vertices  for  searcher  j  E  J .  When  we  describe  an 
arbitrary  searcher,  we  may  omit  the  superscript  and  simply  write  vq  and  V. 

For  any  k  E  T  and  Vi  E  V,  I  =  0,1,2, k,  such  that  {vi-i,vi)  E  £  for 

all  I  =  l,2,...,k,  let  the  sequence  V  =  denote  a  directed  Vo-Vk  subpath.  If 

Vk  E  V,  then  the  directed  VQ-Vk  subpath  is  a  directed  VQ-Vk  path.  In  this  notation, 
a  searcher  flies  from  vq  to  some  Vk  E  V  along  a  directed  vo-Vk  path.  For  a  specific 

searcher  j  E  J,  we  denote  its  directed  v^-Vk  (sub)path  as  Thus  a 

search  plan  for  the  fleet  of  searchers  is  described  asV  =  {V^,V^, 

We  adopt  the  following  target-detection  model.  If  the  target  is  in  cell  c  G 
C  during  time  period  t  E  T  and  only  searcher  j  E  J  is  at  the  same  time  at  a 
waypoint  above  cell  c,  i.e.,  =  c,  then  detection  occurs  with  a  glimpse  detection 

probability  {vl_i,vl,t),  where  vl_i  G  V  is  the  waypoint  of  searcher  j  during  time 
period  t  —  1.  Hence,  the  glimpse  detection  probability  during  time  period  t  depends 
on  the  previous  and  current  waypoints  for  the  searcher.  We  assume  that  the  glimpse 


54 


detection  probability  of  each  searcher  is  mutually  independent.  In  this  notation,  the 
probability  of  detecting  the  target  in  cell  c  G  C  by  any  searcher  during  time  period 
t  G  T  and  no  prior  detections  becomes 

Q(c,t)  =p(c,t)  1-  Yl  {1-  ,  (III.l) 

jej\<l>{vl)=c 

We  refer  to  Q{-,t)  =  [Q{l,t),Q{2,t), . . .  ,Q{C,t)]  as  the  reward  vector  at  time  period 
t. 

Since  p{-,t)  is  the  undetected  target  distribution  at  the  beginning  of  time 
period  t,  it  depends  on  searches  up  to  time  period  t  —  1.  Specifically,  if  p{-,  t)  is  given 
and  each  searcher  searches  a  cell  from  its  waypoint  during  time  period  t  &  T,  the 
undetected  target  distribution  at  the  beginning  of  the  next  time  period  f  +  1  is 

p(.,t  +  l)  =  {p{;t)-Q{;t)]T.  (IIL2) 

Thus,  the  probability  of  detection  for  search  plan  V  is  defined  as 

(ni.3) 

t  =  l  C=1 

We  refer  to  the  problem  of  maximizing  (III. 3)  as  the  multiple-searcher  problem  (MSP). 

B.  ALGORITHMS  FOR  MSP 

We  develop  one  exact  algorithm  and  two  heuristics  to  solve  MSP.  The  exact 
algorithm  is  a  straightforward  extension  of  the  branch-and-bound  (B/B)  algorithm 
in  Chapter  II.  However,  the  B/B  algorithm  for  solving  MSP  requires  an  expanded 
network  structure  to  account  for  multiple  searchers.  Thus,  in  each  time  period  f  G  T, 
we  consider  the  combined  location  Vt  =  (w/,  v/, . . . ,  w/)  of  all  searchers.  We  call 
this  combined  location  a  configuration.  We  treat  a  configuration  in  the  same  way 
we  treated  a  vertex  in  SP,  and  construct  a  time-expanded  (configuration)  network. 
The  B/B  algorithm  for  MSP  is  implemented  in  this  time-expanded  network.  In  that 
network,  a  search  plan  is  simply  a  sequence  of  conhgurations  k  E  T.  We  also 
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use  the  notation  where  ni  =  Thus  we  refer  to  a  search  plan  also  as 

a  “conhguration  path.”  When  the  meaning  is  clear  from  the  context,  we  refer  to  a 
conhguration  path  as  a  path. 

We  present  two  heuristics  to  solve  MSP.  The  hrst  one  is  a  variant  of  the 
expected  detection  heuristic  (HI)  in  [13].  In  HI,  a  conhguration  path  (search  plan) 
is  generated  by  iteratively  extending  the  current  conhguration  subpath  by  the  hrst 
arc  in  the  longest-path,  with  respect  to  the  expected  number  of  detections,  from  the 
current  conhguration  to  the  last  time  period.  We  note  that  a  longest-path  calculation 
is  required  every  time  the  current  conhguration  subpath  is  extended.  Hence,  HI  is 
similar  to  Algorithm  1  with  dynamic  bounds  (see  Chapter  H)  as  they  both  require 
repeated  longest-path  calculations.  As  an  alternative  approach,  we  present  a  similar 
heuristic  corresponding  to  the  static  bound  (see  Chapter  H).  We  denote  our  approach 
the  static  bound  heuristic  (SBH). 

The  second  heuristic  algorithm  is  using  the  Cross-Entropy  (CE)  method  [38, 
39].  The  CE  method  involves  an  iterative  procedure  where  each  iteration  composes 
of  two  phases.  First,  the  method  generates  random  samples  (i.e.,  search  plans  in 
our  case)  according  to  a  probability  distribution.  Second,  the  method  updates  the 
parameters  in  the  probability  distribution  based  on  a  subset  of  the  “best”  samples, 
the  so-called  “elite”  samples.  This  process  increases  the  chances  of  an  optimal  or  near- 
optimal  solution  to  appear  within  the  next  set  of  samples.  This  dissertation  appears 
to  be  the  hrst  one  studying  the  CE  method  in  the  context  of  search  problems. 

We  can  assume  without  loss  of  generality  that  each  searcher’s  path  consists 
of  T  -|-  1  vertices.  For  simplicity  of  notations,  we  also  assume  that:  (1)  there  is  no 
end-point  restriction  for  any  searcher  j  G  JT”,  i.e.,  W  =  V;  (2)  the  glimpse  detection 
probability  g^{v,v',t),  j  E  J'  is  independent  of  v  (i.e.,  g^{v,v',t)  =  g^{v',t));  (3)  the 
airspace  has  only  one  altitude,  i.e.,  there  is  only  one  vertex  corresponding  to  each  cell 
in  the  area  of  interest.  It  is  straightforward  to  generalize  the  proposed  algorithms  to 
account  for  situations  without  these  assumptions. 
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1.  Branch- and-Bound  Algorithm 

In  Chapter  II,  we  presented  a  specialized  branch- and-bound  (B/B)  procednre 
for  solving  SB.  We  ntilize  the  same  enhancement  schemes  to  develop  a  B/B  algorithm 
for  MSP.  However,  the  B /B  algorithm  for  solving  MSP  requires  an  expanded  network 
structure  to  account  for  multiple  searchers.  We  consider  a  directed  graph  (Vc, 
where  Vc  is  the  set  of  possible  conhgurations  and  Ec  is  the  set  of  directed  edges 
representing  possible  transitions  in  one  time  period  between  two  conhgurations.  Thus 
a  search  plan  can  be  described  as  V  =  {VJjiLo  where  (V/-i,Vz)  ^  Ec  for  all  I  = 
1,2,  ...,T.  We  note  that  a  search  plan  is  also  described  (in  the  previous  section)  as 
V  =  {V\  ...Pb  where  =  K}f=o- 

The  description  of  the  B /B  algorithm  for  solving  MSP  follows  the  presentation 
of  the  algorithms  for  SP  presented  in  Chapter  II.  Given  a  conhguration  subpath 
{ViYiZq,  t  G  T,  let  /C(t)  be  the  set  of  triplets  of  the  form  q{Vt,t))  representing 
extensions  of  {V;};=o  fo  be  explored.  The  hrst  element  V)  refers  to  the  next 
conhguration  to  form,  the  second  element  t  is  the  time  period  to  form  the  conhguration 
Vt,  and  the  third  element  q(yt,t)  is  an  upper  bound  on  the  probability  of  detection 
along  any  conhguration  path  that  starts  with  the  conhguration  subpath  {V)}f=o-  The 
upper  bound  g(V),t)  consists  of  three  parts.  Let  be  an  upper  bound  on 

the  probability  of  detection  during  time  periods  t  +  l,t  +  2,  ...,T,  given  that  the 
conhguration  of  time  t  is  Vt,  and  no  detection  occurs  along  the  conhguration  subpath 
The  two  other  parts  are  the  probability  of  detection  on  the  conhguration 
subpath  {Vi}\Zq  and  the  probability  of  detection  during  t.  Hence, 

q{Vu  t)  =  qiiVYzY  +  J2Qivt)  +  d,( V„  t).  (HI.4) 

cec 

We  also  let  q  denote  the  largest  detection  probability  found  so  far  among  all  the  ex¬ 
amined  conhguration  paths. 
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Consider  a  configuration  subpath  {V;}*^q,  t  eT,  and  let  Pg{-,t)  be  the  unde¬ 
tected  target  distribution  after  search  along  i.e., 

Pgi-,^)  =  (III.5) 

For  any  integer  s  >  t,  s  E  T,  we  also  define 

Pr{;s;t)=pg{;t)r-K  (III.6) 

As  seen,  pr{c,s;t)  is  the  probability  that  the  target  is  in  cell  c  at  time  period  s  >  t 
and  there  was  no  detection  during  search  along  the  subpath  Hence,  target 

distribution  pr{-,s]t)  at  time  period  s  >  t  ignores  the  effect  of  search  after  time 
period  t.  If  the  subpath  is  {Vq},  i.e.,  t  =  0,  we  dehne  for  notational  convenience 

pr(-,s;0)=p(-,l)r-\  (III.7) 

for  any  s  >  0,  s  G  T.  Moreover,  we  dehne  pr{c,  t;t)  =  0  for  all  c  G  C  and  t  =  0, 1, ...,  T. 

We  construct  a  time-expanded  conhguration  graph  (A/'c,  Ac)  from  the  graph 
(Vc,  £c)  in  the  same  manner  as  the  development  of  the  time-expanded  network  for 
SP  in  Chapter  II  (for  details,  we  refer  to  Section  B  in  Chapter  II).  For  any  integer 
k  <  T  +1  and  nodes  ni  =  {Vi,l)  E  Me,  I  =  0,1,. ..,k,  such  that  (nj_i,nj)  G  Ac  for 
all  I  =  l,2,...,k,  we  let  the  sequence  {ni}f^Q  denote  a  conhguration  subpath  in  the 
time-expanded  conhguration  graph  {Me,  Ac). 

For  some  t  E  {0,1,...,T  —  1},  suppose  that  a  conhguration  subpath 
in  the  graph  (Vc,  Sc)  is  given.  We  endow  each  arc  (n,  n')  =  ((V,  s),  {V ,  s  -f- 1))  G  Ac, 
s  =  t,t  +  l,...,T  —  1,  in  the  time-expanded  conhguration  graph  {Me,  Ac)  with  a 
“reward” 


- 

- 

Pr{c’,s  +  l;t)  - 

QMs)T{c,c') 

n  -  9^ s  +  ^) 

c'GC 

- 

ce7^(c') 

(IIL8) 


where  F(c,  c')  is  the  c-c'  element  of  the  Markov  transition  matrix  F.  Thus  this  time- 
expanded  conhguration  network  {Me,  Ac)  has  the  same  structure  as  the  time-expanded 
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network  (A/*,  A)  of  Chapter  II  and  the  longest  paths  in  this  network  provides  the  static 
bound  dolVt,  t).  The  B/B  for  MSP  is  implemented  in  this  time-expanded  conhguration 
network  {Nc,Ac). 

We  consider  generalizations  of  directional  static  bound  and  network  reduction 
from  Chapter  II.  To  obtain  a  directional  static  bounds,  we  would  need  to  generate  a 
node-and-time  expanded  “conhguration”  network  corresponding  to  the  node-and-time 
expanded  network  for  SP.  Since  the  number  of  arcs  in  the  node-and-time  expanded 
conhguration  network  becomes  extremely  large  (e.g.,  a  problem  instance  with  5  by  5 
cells,  3  searchers,  and  time  horizon  7  has  more  than  one  hundred  million  arcs),  we 
realize  that  the  directional  static  bound  technique  is  not  computationally  practical  for 
MSP.  Hence,  we  adopt  the  static  bound  in  our  B/B  algorithm  for  MSP.  The  difficulty 
of  network  size  occurs  also  in  the  time-expanded  conhguration  network  when  we 
consider  large-size  problems,  which  means  our  B/B  algorithm  is  practical  only  for 
small-size  problems. 

It  is  straightforward  to  generalize  the  network-reduction  technique  based  on 
the  initial  positions  of  the  searchers  and  the  target  of  Chapter  II  to  MSP.  The  required 
“conhgurations”-of-hrst-contact  (14,  s)  G  A//  are  now  dehned  as  the  hrst  conhgura¬ 
tion  where  at  least  one  searcher  can  obtain  contact  with  the  target.  Thus  our  B/B 
algorithm  utilizes  this  network  reduction  technique.  The  resulting  algorithm  takes 
the  following  form: 

Algorithm  4  (Solves  MSP). 

Input:  A  time-expanded  conhguration  network  (A//,  Ac)  with  nodes  n  G  A/” 
and  arcs  (n, n')  G  A  where  n  =  {V,t  —  1),  n'  =  Arc  lengths  c„^„/ 

in  (A//,  Ac),  the  initial  target  distribution  p{-,  1),  Markov  transition  matrix 
P,  and  upper  bound  q. 

Output:  An  optimal  conhguration  path  V*  =  {V)}^o  value  q* . 

Step  0.  Calculate  the  static  bound  do(y,t)  for  all  t  G  T  and  V  G  Vc,  and 
calculate  a  lower  bound  q.  Set  i  =  0,/C(i)  =  {(Vq;  0,  cx))},  where  Vq  is 
the  initial  conhguration.  Implement  network  reduction  procedure  R5  (see 
below) . 
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Step  1.  If  /C(f)  is  empty,  replace  ihy  i  —  1.  Else,  go  to  Step  3. 

Step  2.  \i  i  =  0,  stop:  the  last  saved  configuration-path  is  optimal  and  q  is 
its  probability  of  detection.  Else,  go  to  Step  1. 

Step  3.  Remove  from  lC{i)  the  triplet  (yt,t,q(yt,t))  with  the  largest  bound 

Step  4.  If  qiVtjt)  <  q,  go  to  Step  1.  (Current  subpath  is  fathomed.) 

Step  5.  If  f  =  0,  replace  ihy  i  +  1,  and  go  to  Step  3.  (/C(l)  is  populated  in 
the  network  reduction  procedure.) 

Step  6.  If  t  R,  then  for  each  conhgnration  R  G  .^(V)),  calculate  q(V^t  1) 
from  (III. 4)  using  bounds  do{V,t  +  1)  and  add  {V,t  -|-  l,q{V,t  +  1))  to 
)C{i  +  1).  Replace  f  by  f  -|-  1,  and  go  to  Step  3.  Else,  let  q  =  qiVt^t)  and 
save  the  incumbent  conhgnration  path  {V/}^g,  and  go  to  Step  1. 

Network  Reduction  Procedure  R5. 

Input:  A  time-expanded  conhgnration  network  (A/),, .Ac),  the  initial  target 
distribution  p(-,  1),  Markov  transition  matrix  E,  and  upper  bound  q. 

Output:  The  conhgurations-of-hrst-contacts  (K,<s)  which  are  not  dominated 
by  others. 

Step  1.  Find  all  conhgurations-of-hrst-contact  (IA,s)  G  Me-  If  none  exist 
with  s  >  1,  then  stop. 

Step  2.  For  each  (14,  s),  calculate  g(I4,  s)  =  I]cec  Q(c,  -s)  -|-  (io(f4,  -s), 

Step  3.  Eliminate  all  conhgurations-of-hrst-contact  (14,  s)  with  g(I4,  s)  <  q. 

Step  4.  For  each  conhgurations-of-hrst-contact  (14,  s)  not  eliminated,  store 
the  triplet  (14,  s,  qiVs,  s))  in  /C(l). 

We  implement  Algorithm  4  and  examine  its  performance  on  small-scale  in¬ 
stances  of  MSP.  The  problem  instances  are  as  follows:  An  area  of  interest  (AOI) 
consists  of  5  by  5  cells.  The  number  of  searchers  is  1,  2,  and  3.  The  searchers,  which 
have  a  constant  glimpse  detection  probability  0.6,  operate  in  the  AOI  and  their  moves 
are  restricted  to  vertically  or  horizontally  adjacent  cells.  The  target  remains  in  the 
current  cell  with  probability  0.6  or  moves  to  one  of  the  vertically  or  horizontally  ad¬ 
jacent  cells  with  probability  0.4.  The  diherent  moves  are  equally  likely.  All  searchers 
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depart  the  upper-left  corner  cell,  and  the  target  starts  at  the  center  of  the  AOI.  The 
time  horizon  is  from  5  to  10.  We  carry  out  all  experiments  on  the  same  computation 
platform  as  in  Chapter  II. 

Table  8  reports  the  run  times  for  different  combinations  of  number  of  searchers 
and  time  horizons.  Column  2,  4  and  6  of  Table  8  show  the  run  times  of  Algorithm 
4  for  the  case  of  1,  2  and  3  searchers,  respectively.  For  the  case  with  3  searchers, 
the  run  time  increases  exponentially  with  increasing  time  horizon,  and  the  case  with 
time  horizon  10  cannot  be  solved  in  several  days.  The  case  with  3  searchers  and 
the  time  horizon  9  takes  about  5  days  to  guarantee  an  optimal  solution,  however, 
the  optimal  solution  is  found  after  only  174.78  seconds.  The  remaining  time  is  spent 
proving  optimality  by  fathoming  other  possible  conhguration  paths.  Table  9  presents 
the  optimal  probability  of  detection  obtained  by  Algorithm  4.  Column  2,  5  and  8  of 
Table  9  show  the  optimal  probability  of  detection  for  the  case  of  1,  2  and  3  searchers, 
respectively.  Figure  9  illustrate  an  optimal  search  plan  for  the  case  with  3  searchers 
and  time  horizon  9. 


Time 

Horizon 

1  searcher 
B/B  SBH 
(sec.)  (sec.) 

2  searchers 
B/B  SBH 
(sec.)  (sec.) 

3  searchers 

B/B  SBH 

(sec.)  (sec.) 

5 

0.00 

0.00 

0.02 

0.00 

4.17 

1.84 

6 

0.00 

0.00 

0.03 

0.02 

5.67 

2.17 

7 

0.00 

0.00 

0.06 

0.02 

58.89 

2.72 

8 

0.00 

0.00 

0.44 

0.02 

4,341.51 

2.91 

9 

0.00 

0.00 

5.19 

0.03 

426,874.92 

3.25 

10 

0.00 

0.00 

75.77 

0.03 

N/A 

3.70 

Table  8.  Run  times  for  the  branch-and-bound  algorithm  (B/B)  and  the  static  bound 
heuristic  (SBH)  on  5  by  5  cell  search  problem  with  1-3  searchers  and  time  horizon 


T  =  5-10. 


In  view  of  Tables  8  and  9,  it  is  clear  that  Algorithm  4  is  practical  only  for 
small-size  problems  (with  few  searchers,  short  time  horizon,  and/or  small  area  of 
interest).  As  noted  earlier,  the  time-expanded  configuration  network  (necessary  for 
implementing  the  B/B  procedure)  cannot  be  generated  for  large-size  problems  since 
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Time 

Horizon 

1  searcher 
B/B  SBH 

Rel. 

gap 

2  searchers 
B/B  SBH 

Rel. 

gap 

3  searchers 
B/B  SBH 

Rel. 

gap 

5 

0.306483 

0.306483 

0.00 

0.474213 

0.474213 

0.00 

0.579710 

0.579710 

0.00 

6 

0.351647 

0.351241 

0.00 

0.535954 

0.521669 

0.03 

0.643001 

0.622074 

0.03 

7 

0.389043 

0.380220 

0.02 

0.581175 

0.561550 

0.03 

0.691865 

0.679234 

0.02 

8 

0.416987 

0.404325 

0.03 

0.618416 

0.574542 

0.07 

0.728375 

0.711876 

0.02 

9 

0.444506 

0.426829 

0.04 

0.647400 

0.620582 

0.04 

0.754400 

0.739376 

0.02 

10 

0.465594 

0.438671 

0.06 

0.673168 

0.648007 

0.04 

N/A 

0.762183 

N/A 

Table  9.  Probability  of  detection  obtained  by  the  branch- and-bound  algorithm  (B/B) 
and  the  static  bound  heuristic  (SBH)  on  5  by  5  cell  search  problem  with  1-3  searchers 
and  time  horizons  T  =  5-10,  and  relative  gap  between  SBH  and  B/B  solutions. 


the  extremely  large  number  of  data  (arcs)  can  not  be  stored.  Thus,  to  solve  large-size 
problems,  alternative  heuristic  approaches  may  be  necessary. 

2.  Static  Bound  Heuristic 

Dell  et  ah  [13]  present  seven  algorithms  (one  exact  procedure  and  six  heuris¬ 
tics)  to  solve  MSP.  Among  their  heuristics,  the  expected  detection  heuristic  (HI) 
performs  best  under  a  broad  range  of  conditions.  As  described  in  the  beginning  of 
this  section,  HI  constructs  a  search  plan  (conhguration  path)  by  extending  a  config¬ 
uration  subpath  one  arc  at  the  time.  Each  extension  is  determined  by  a  longest-path 
calculation  on  a  time-expanded  conhguration  network  where  arc  lengths  are  given  by 
the  expected  number  of  detections. 

We  present  an  alternative  approach,  the  static  bound  heuristic  (SBH),  which  is 
similar  to  HI  in  [13]  but  requires  only  one  longest-path  calculation  and  uses  improved 
arc  lengths  as  described  below.  In  Algorithm  4  (i.e.,  the  B/B  procedure  for  MSP),  we 
calculate  the  static  bounds  in  advance  of  the  main  calculations.  The  bound  calculation 
corresponds  to  a  longest  path  in  the  time-expanded  conhguration  network  (A//,  Me) 
with  arc  length  given  by  (HI. 8).  That  longest  path,  which  specihes  a  search  plan, 
is  the  solution  provided  by  SBH.  We  observe  that  the  longest  path  calculations  in 
HI  amounts  to  hnding  a  conhguration  path  with  the  largest  expected  number  of 


62 


— 

i 

k — 

— ► 

r 

> 

Searchers 


Figure  9.  An  optimal  search  plan  for  the  5  by  5  cell  search  problem  with  3  searchers 
and  the  time  horizon  9. 

detections.  In  contrast,  the  longest  path  calculation  in  SBH  uses  arc  lengths  (III. 8) 
which  effectively  amounts  to  hnding  the  conhguration  path  with  the  largest  expected 
number  of  detections  while  at  each  node  along  the  path  accounting  for  the  effect  of 
search  at  the  previous  node.  As  demonstrated  in  [32]  and  discussed  in  Section  1  of 
Chapter  II,  accounting  for  the  previous  node  signihcantly  improves  the  conhguration 
path  found  by  the  longest  path  calculation. 

Column  3,  5  and  7  of  Table  8  report  the  run  times  for  the  cases  of  1,  2  and 
3  searchers,  respectively.  Since  SBH  corresponds  to  a  longest-path  problem  in  an 
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acyclic  network,  the  run  times  tend  to  be  extremely  short.  Column  3,  6  and  9  of 
Table  9  show  the  obtained  probability  of  detection,  and  column  4,  7  and  10  describe 
the  relative  gap  to  the  corresponding  optimal  probability  of  detection.  For  each  case, 
the  obtained  detection  probability  is  quite  close  to  the  corresponding  optimal  value. 
Thus  SBH  performs  well  for  these  problem  instances.  However  the  performance  tends 
to  be  worse  as  the  problem  size  becomes  large.  Moreover,  this  SBH  is  not  available 
for  large-size  problem  since  the  time-expanded  conhguration  cannot  be  generated. 

3.  Cross-Entropy  Method 

The  cross-entropy  (CE)  method  was  developed  by  Rubinstein  [38,  39]  for  solv¬ 
ing  rare  event  simulations  and  combinatorial  optimization  problems.  The  CE  method 
derives  its  name  from  the  cross-entropy  (or  Kullback-Leibler)  distance  between  two 
probability  distributions  (e.g,  an  optimal  importance  sampling  distribution  and  an 
estimated  distribution).  The  CE  method  is  an  iterative  procedure  consisting  of  two 
steps:  (1)  random  samples  (i.e.,  search  plans  for  our  case)  are  generated  according  to 
a  parameterized  probability  distribution,  and  (2)  the  generated  search  plans  are  eval¬ 
uated  using  the  objective  function  (i.e.,  detection  probability),  and  the  parameters 
of  the  sampling  distribution  are  updated  based  on  the  ’’elite”  samples  in  a  manner 
which  increases  the  possibility  that  an  optimal  or  near-optimal  solution  appears  in 
the  next  iteration. 

The  CE  method  is  a  global  search  heuristic,  and  is  somewhat  similar  to  ge¬ 
netic  algorithms.  However  the  CE  method  has  a  simpler  scheme  for  the  change  of 
population  generation  parameter  compared  to  genetic  algorithms.  Boer  et  ah  [36] 
compare  in  details  the  CE  method  and  other  heuristics  such  as  simulated  annealing, 
genetic  algorithm,  and  ant  colony  method.  The  CE  method  has  been  applied  to  a 
large  number  of  combinatorial  optimization  problems  [39,  40,  33,  44,  36,  12]  but  this 
dissertation  appears  to  be  the  first  one  applying  the  CE  method  to  search  problems. 
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a.  CE  Algorithms 

In  the  multiple-search  problem  (MSP),  a  possible  search  plan  is  de¬ 
scribes  as  a  sequence  of  conhgurations  V  =  {Vz}^q  where  (V)_i,V/)  G  £c  for  all 
I  =  Thus,  a  crude  way  to  hud  a  search  plan  is  as  follows.  Let  EiV)  = 

{V  G  Vc|  (Id,  V')  G  £c}  be  the  forward  star  of  conhguration  V.  Start  at  the  given  ini¬ 
tial  conhguration  Vq  =  (^0,^0;  •  •  •  Select  an  arbitrary  conhguration  from  JP(Vo) 

and  denote  it  Vi.  Next,  choose  an  arbitrary  conhguration  from  EiVi)  and  denote  it 
V2.  The  same  process  is  repeated  until  the  conhguration  Vt  is  obtained. 

The  CE  method  takes  similar  steps,  but  selects  the  next  conhguration 
according  to  a  probability  distribution.  For  each  node  n  G  Me  in  the  time-expanded 
conhguration  network  (A/);, Me),  we  dehne  a  probability  distribution  an,n'  over  the 
outgoing  arcs  (n,  n')  G  Ac  (i.e.,  J^(n,n')GAc  Mn'  =  !)•  Then,  given  a  conhguration 
subpath  k  E  T,  the  conhguration  subpath  is  extended  by  randomly  selecting 

an  arc  {nk,n')  G  Ac  according  to  the  probability  distribution  <Jnk,n'- 

Clearly,  if  an,n'  is  degenerate  at  each  node  n  G  M^  i.e.,  there  exists  an  n' 
such  that  (n,  n')  G  Ac  and  an,n'  =  1,  then  the  probability  distribution  uniquely  speci- 
hes  a  conhguration  path.  The  CE  methods  aims  to  iteratively  update  the  probability 
distribution  so  that  it  converges  to  a  degenerate  distribution  specifying  the  optimal 
conhguration  path  of  MSP.  We  formalize  this  approach  in  the  next  algorithm. 


Algorithm  5  (Basic  CE  Algorithm). 

Parameters.  Sample  size  M,  elite  size  stopping  parameter  s  (e.g., 

s  =  5),  and  smoothing  parameter  /3. 


Step  0.  Calculate  static  bounds  do(n)  for  all  n  E  Me  and  a  lower  bound  q. 
Set  =  q  and  i  =  1. 


Step  1.  For  all  n  E  Me  and  n'  E  A’(n),  dehne  probability  distribution 
such  that 

(i)  _  do(n') 


^£2in'(iT(n)  df)(n  ) 


(III.9) 


Step  2.  Generate  M  search  plans  based  on  probability  distribution  a 


(i) 

n,n'  * 
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Step  3.  Choose  elite  samples  among  the  generated  search  plans,  and 

calculate  the  detection  probability  of  the  best  elite  If  >  g,  then 
g  =  g*^*) 

Step  4.  If  gW  =  g(*-i)  =  . . .  =  gh“^),  then  stop.  Else  go  to  Step  5. 

Step  5.  Update  probability  distribution  using  the  elite  samples 

and  parameter  f3  as  described  below.  Replace  i  by  i  +  1  and  go  to  Step  2. 

In  Step  0,  the  static  bounds  d^in)  are  calculated  for  all  node  n  G  A/'c  in 
the  time-expanded  conhguration  network.  At  this  time,  we  obtain  a  lower  bound  g 
of  the  optimal  detection  probability  (i.e.,  the  detection  probability  along  the  conhg¬ 
uration  path  corresponding  to  do(no)).  The  best  detection  probability  initially  (0th 
iteration),  denoted  g*^°\  is  simply  g.  In  Step  1,  we  dehne  the  initial  probability  dis¬ 
tribution  using  the  static  bounds.  Next,  in  Step  2,  based  on  the  probability 

distribution  we  generate  M  search  plans  Pi,  7^2,  •  •  • ,  Pm  randomly.  After  that, 
in  Step  3,  the  M  plans  are  evaluated  using  the  objective  function  (i.e.,  detection  prob¬ 
ability),  and  the  best  performing  elite  samples  are  extracted.  If  the  detection 

probability  of  the  current  best  elite  sample  g*^^^  is  better  than  the  current  lower  bound 
g  (the  best  detection  probability  so  far),  g  is  updated.  Step  4  stops  the  algorithm  if 
the  best  values  found  in  the  previous  several  iterations  were  identical.  We  hope  that 
the  repetition  of  values  is  due  to  a  near  degenerate  probability  distributions  and 
that  the  best  found  search  plan  is  close  to  the  optimal  one.  However,  that  is  not  guar¬ 
anteed.  If  the  stopping  criterion  is  not  satished  in  Step  4,  probability  distributions 
are  updated  using  the  current  elite  samples  as  follows.  Let  be  the  number 

of  times  one  of  the  elite  samples  use  the  arc  (n,  n')  in  the  current  iteration.  The 
probability  distribution  is  updated  using  both  the  current  probability  distributions 
and  contribution  of  the  elite  samples  by  setting 

= fi  +  (1  -  -  (iii.io) 

^n'£j^{n)  ^n,n' 

where  0.4  <  P  <  0.9  as  suggested  in  [36]  based  on  empirical  studies. 


66 


We  note  that  the  parameters  (M,  s,  j3)  must  be  specihed  to 

implement  Algorithm  5.  Boer  et  ah  [36]  suggest  that  the  sample  size  M  is  chosen 
according  to  the  size  of  the  problem  (e.g.,  M  =  rm,  where  m  is  the  number  of  arcs  in 
the  time-expanded  conhguration  network,  and  r  is  a  constant),  and  the  elite  sample 
size  is  taken  to  be  approximately  O.OIM.  Boer  et  ah  [36]  also  report  that  an  adaptive 
choice  of  the  parameter  settings  speeds  up  the  convergence.  Thus,  we  implement  nu¬ 
merical  tests  for  a  variety  of  cases  and  develop  an  adaptive  CE  algorithm.  Specifically, 
we  set  the  elite  sample  size  and  smoothing  parameter  as  =  O.OIM  and  (3  =  0.4, 

respectively.  Furthermore,  we  vary  the  sample  size  M  in  the  range  between 
and  adaptively.  M™'*"'  is  set  according  to  the  size  of  the  problem  instance,  and 

j^max  _ 

The  CE  method  has  asymptotic  convergence  properties.  Specifically, 
under  the  assumption  that  the  optimal  solution  is  unique,  Costa  et  ah  [3]  give  neces¬ 
sary  and  sufficient  conditions  under  which  the  optimal  solution  is  approached,  with 
probability  one,  as  the  number  of  iterations  goes  to  inhnity.  For  MSP,  the  CE  method 
provides  good  solutions  (as  shown  below  in  numerical  tests).  However  it  does  not 
guarantee  an  optimal  solution  since  MSP  often  has  multiple  optimal  solutions  (e.g.,  if 
two  searchers  are  identical)  and  since  the  CE  method  is  implemented  with  a  heuristic 
stopping  criterion.  The  multi-extremal  property  of  the  MSP  solutions  may  provide 
the  following  undesirable  condition.  In  each  iteration,  probability  distributions  are 
updated  using  the  elite  samples  including  the  best  elite.  Thus,  in  the  next  iteration, 
at  least,  the  best  elite  of  the  previous  iteration  is  expected  to  appear  in  the  samples 
if  the  sample  size  is  large.  However,  the  multi-extremal  property  may  bring  unstable 
probability  distributions  and  the  current  best  elite  may  be  worse  than  that  in  the 
previous  iteration.  In  this  case.  Algorithm  5  is  unlikely  to  stop  in  reasonable  time. 
For  that  undesirable  situation,  we  intentionally  stop  the  algorithm  and  provide  the 
best  solution,  but  label  it  as  “unreliable.” 
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In  view  of  the  above  discussion,  we  develop  a  specialized  adaptive  CE 
algorithm  to  solve  the  MSP.  In  this  algorithm,  we  refer  to  the  sample  size  at 
iteration  as  Mi. 


Algorithm  6  (Adaptive  CE  Algorithm). 

Parameters.  Minimum  and  maximum  sample  sizes  and  and 

smoothing  parameter  13. 


Step  0.  Calculate  static  bounds  do{n)  for  all  n  G  Me  and  lower  bound  q.  Set 
=  q  and  i  =  1. 


Step  1.  For  all  n  G  Me  and  n'  G  JP(n),  define  probability  distribution  a. 


(i) 

n,n' 


such  that 


= 

^n,n' 


doin') 


(IIMl) 


^2in' (iT(n)  d^ijl  ) 

Set  Mi  =  M^^^. 

Step  2.  Generate  Mi  search  plans  based  on  probability  distribution  u^nn'- 


Step  3.  Choose  =  O.OlMj  elite  samples  among  the  generated  search 

plans,  and  calculate  the  detection  probability  of  the  best  elite  If  > 
g,  then  g  =  g^*^ 

Step  4.  If  gh)  >  gh-P^  update  using  the  elite  samples.  Replace 

ihy  i  +  1.  Set  Mi  =  M™*”  and  go  to  Step  2.  Else  go  to  Step  5. 


Step  5.  If  the  best  detection  probabilities  are  identical  for  three  iterations 
(i.e.,  g^*^  =  gh“^)  =  then  stop.  (The  obtained  solution  is  considered 

reliable).  Else  go  to  Step  6. 

Step  6.  If  Mj  <  increase  the  sample  size  Mi  =  Mi  + 

and  go  to  Step  2.  Else  go  to  Step  7. 

Step  7.  If  Mk  =  Mi_i  =  ...  =  Mi_^,  then  stop.  (The  obtained  solution 
is  considered  unreliable).  Else  update  using  the  latest  elite 

samples.  Replace  ihy  i  +  1.  Mj  =  M™'*"'  and  go  to  Step  2 


Algorithm  6  is  based  on  the  time-expanded  conhguration  network, 
which  is  problematic  for  large-scale  problems.  However,  we  overcome  this  difficulty 
by  considering  a  time-expanded  network  for  each  searcher.  For  each  searcher  j  G  77, 
we  construct  a  time-expanded  network  as  in  Chapter  II.  For  each  Q\  we  dehne  the 
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probability  distributions  We  then  generate  a  search  plan  by  generating  a  path 

=  {n/}^Q  for  each  searcher  j  according  to  the  probability  distributions  The 
search  plan  is  simply  V  =  We  update  the  probability  distributions 

j  G  J',  by  considering  the  joint  search  effect  of  all  searchers.  Specihcally,  we 
determine  as  before  the  elite  search  plan  and  dehne  be  the  number  of  times 

one  of  the  elite  samples  use  the  arc  (n,  n')  for  search  j  in  the  current  iteration. 
The  probability  distribution  for  each  searcher  j  is  then  updated  by  setting 


_  a 

^n,n'  —  P 


^n,n' 


/i(b 

2^n'£T(n)  ^n,n' 


+  (1  ~  /^)t 


J'h) 

n,n'  ‘ 


(III.12) 


An  advantage  of  this  approach  is  that  we  can  treat  large-scale  problems  since  we 
avoid  generating  the  time-expanded  conhguration  network.  We  summarize  this  “de¬ 
composition”  approach  next. 


Algorithm  7  (Adaptive  CE  Algorithm  with  Decomposition). 

Parameters.  Minimum  and  maximum  sample  sizes  and  and 

smoothing  parameter  /?. 

Step  0.  Generate  the  time-expanded  network  ,  j  G  J .  In  each  calculate 
static  bounds  d^iji)  for  all  nodes  n  G  .  Generate  a  search  plan  V  by 
collecting  each  searcher’s  path  corresponding  to  the  static  bound  dh^no). 
Galculate  detection  probability  g('P)  from  (III. 3).  Set  q  =  =  g('P)  and 

Z  =  1. 

Step  1.  For  j  G  JT”,  define  the  probability  distributions  by 


Set  Mi  = 


^  doin') 
En'er(n)di{n^y 


(III.13) 


Step  2.  Gonstruct  Mi  search  plans  by  generating  each  searcher’s  path  based 
on  probability  distribution  For  each  search  plan,  V  calculate  the 

detection  probability  qiV)  from  (III. 3). 

Step  3.  Ghoose  =  O.OlMj  elite  search  plans  among  the  generated 

search  plans,  and  set  the  detection  probability  of  the  best  elite  as  qd\ 
If  qd')  >  g,  then  q  =  qd) 
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Step  4.  If  >  update  Vj  G  J  using  searcher  j’s  path  con¬ 
tributing  to  elite  samples.  Replace  f  by  i  -|-  1.  Mj  =  M™'*"'  and  go 

to  Step  2.  Else  go  to  Step  5. 

Step  5.  If  the  best  detection  probabilities  are  same  for  three  iterations  (i.e., 
gW  =  g(*-i)  =  then  stop.  (The  obtained  solution  is  considered 

reliable).  Else  go  to  Step  6. 

Step  6.  If  Mj  <  increase  the  sample  size 

and  go  to  Step  2.  Else  go  to  Step  7. 

Step  7.  If  Mfc  =  Mj_i  =  . . .  =  Mj_5,  then  stop.  (The  obtained  solution  is 
considered  unreliable).  Else  update  j  G  77,  using  (III.  12).  Replace 

f  by  f  -f  1.  Set  Mi  =  M™*”  and  go  to  Step  2 

In  Algorithm  7,  for  each  j  G  77,  the  probability  distributions 
are  rather  poor  in  the  initial  iterations  compared  to  the  situation  in  Algorithm  6. 
However,  the  distributions  are  gradually  improved  as  the  cooperative  search  effort  is 
accounted  for  when  determining  the  elite  search  plans. 

The  decomposition  approach  is  practical  for  large-size  problems  since  it 
only  requires  the  time-expanded  network  for  each  single  searcher.  Thus,  the  approach 
is  easily  be  implemented  using  the  node-and-time  expanded  network  with  the  poten¬ 
tial  for  more  effective  generation  of  paths  (see  Chapter  II  for  a  discussion  about  the 
advantage  of  the  note-and-time  expanded  network  over  the  time-expanded  network). 
The  size  of  the  node-and-time  expanded  network  is  larger  than  time-expanded  net¬ 
work,  but  it  is  substantially  smaller  than  the  time-expanded  conhguration  network. 
Since  the  algorithm  using  the  node-and-time  expanded  network  is  essentially  identical 
to  Algorithm  7  (except  for  using  the  node-and-time  expanded  network),  we  omit  to 
explicitly  describe  that  algorithm  but  refer  to  it  as  Algorithm  7N. 

h.  Numerical  Tests 

We  implement  Algorithm  6,  Algorithm  7,  and  Algorithm  7N  using  the 
twelve  problem  instances  listed  in  Table  10.  The  problem  instances  are  similar  to  the 
ones  in  Tables  8  and  9,  and  they  have  the  same  assnmptions  and  parameter  settings 
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(e.g.,  start  positions,  glimpse  detection  probability  and  target  stationary  probability, 
etc).  The  minimum  sample  size  M™*"  is  decided  according  to  the  size  of  problem 
instance,  see  columns  5  and  6  of  Table  10  where  m  is  the  number  of  arcs  in  the 
network  corresponding  to  the  algorithms.  As  before,  we  set  =  5M™*"  and 

(3  =  0.4.  Since  Algorithm  6  uses  the  time-expanded  conhguration  network,  it  is 
available  only  for  small  problem  instances,  i.e.,  cases  A2,  A3,  B2,  and  B3.  The  CE 
heuristics  (Algorithms  6,  7,  and  7N)  are  compared  with  Algorithm  4  (i.e.,  the  B/B 
algorithm).  When  available,  we  compare  the  optimal  value  from  Algorithm  4  with  the 
search  plan  obtained  by  the  CE  heuristics.  We  also  compare  with  incumbent  solution 
available  in  Algorithm  4  at  the  point  in  time  when  the  CE  heuristics  terminate. 


Case 

Area 

size 

Time 

horizon 

Number  of 
searchers 

Min.  sample  size 
Algo.  6  Algo.  7/7N 

A2 

5  by  5 

9 

2 

O.lm 

10000m 

A3 

5  by  5 

9 

3 

0.01m 

10000m 

B2 

5  by  5 

18 

2 

O.lm 

10000m 

B3 

5  by  5 

18 

3 

0.01m 

10000m 

C2 

15  by  15 

18 

2 

N/A 

1000m 

C3 

15  by  15 

18 

3 

N/A 

1000m 

D1 

15  by  15 

27 

1 

N/A 

100m 

D2 

15  by  15 

27 

2 

N/A 

100m 

D3 

15  by  15 

27 

3 

N/A 

100m 

El 

15  by  15 

28 

1 

N/A 

100m 

FI 

15  by  15 

29 

1 

N/A 

100m 

G1 

15  by  15 

30 

1 

N/A 

100m 

Table  10.  Test  cases  with  Algorithm  6,  Algorithm  7  and  Algorithm  7N  with  area  size, 
time  horizon,  number  of  searchers,  and  minimum  sample  size,  m  is  the  number  of 
arcs  in  the  network  used  in  the  corresponding  algorithm. 
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(a)  Single  searcher 

Tables  11  and  12  report  the  numerical  results  for  problem  instances  with 
a  single  searcher.  Columns  2  and  3  in  Table  11  present  the  detection  probabilities 
obtained  by  the  CE  heuristics  and  column  4  shows  the  corresponding  optimal  solution 
obtained  by  Algorithm  4.  We  hnd  that  Algorithms  7  and  7N  have  comparable  solution 
quality.  The  last  column  in  Table  11  reports  the  probability  of  detection  of  the 
incumbent  search  plan  of  Algorithm  4  at  the  time  Algorithm  7  terminates.  For 
example,  for  case  Dl,  the  best  detection  probability  available  to  Algorithm  4  at  time 
50.98  seconds  (see  Table  12)  is  0.304046.  Table  12  reports  the  corresponding  run 
times  (columns  2  and  3  for  Algorithms  7  and  7N,  respectively,  and  column  5  for 
Algorithm  4).  In  Table  12,  we  have  also  included  the  time  at  which  Algorithm  4  Ends 
the  optimal  solution  (but  not  yet  proven  optimal),  see  column  4.  The  run  time  of 
Algorithm  4  is  clearly  slower  than  that  of  the  CE  heuristics.  We  also  observe  that  the 
run  time  of  Algorithm  7  is  much  faster  than  that  of  Algorithm  7N.  In  a  comparison 
of  columns  2,  3,  and  5  in  Table  11,  we  find  that  the  CE  heuristics  terminates  with 
a  better  search  plan  than  what  is  available  from  Algorithm  4  at  the  same  time.  The 
solution  quality  of  the  CE  heuristics  is  quite  good  for  all  the  instances  examined. 


Case 

CE  heuristic 

Algo.  7  Algo.  7N 

Branch-and-bound  (Algo.  4) 
Optimal  Early  term. 

Dl 

0.305254 

0.305254 

0.305254 

0.304046 

El 

0.311225 

0.313101 

0.313101 

0.307411 

FI 

0.320277 

0.320716 

0.320719 

0.313043 

G1 

0.325968 

0.325968 

0.327823 

0.320131 

Table  11.  Probability  of  detection  for  problem  instances  with  a  single  searcher. 


(h)  Two  searchers 

Tables  13  and  14  describe  the  computational  results  for  problem  in¬ 
stances  with  two  searchers.  The  columns  show  the  same  content  as  the  corresponding 
columns  in  Tables  11  and  12.  Algorithm  4  performs  well  for  the  smallest  problem 
instance  (A2),  but  is  not  available  for  large  instances.  Algorithms  7  and  7N  obtain 
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CE  heuristic 

Branch-and-bound  (Algo.  4) 

Algo.  7 

Algo.  7N 

Found  optimal 

Proved  optimal 

Case 

(sec.) 

(sec.) 

(sec.) 

(sec.) 

D1 

50.98 

226.08 

151.14 

752.32 

El 

47.28 

229.17 

536.40 

2732.36 

FI 

79.26 

302.88 

1955.64 

9690.79 

G1 

60.17 

309.08 

6255.47 

34820.26 

Table  12.  Run  times  on  problem  instances  with  a  single  searcher. 


solutions  in  practical  times  and  their  solution  quality  is  better  than  that  of  Algorithm 
6.  In  addition,  they  can  obtain  solutions  for  large  problems.  Hence,  it  appears  that 
Algorithms  7  and  7N  are  superior  to  Algorithm  6  for  problem  instances  with  two 
searchers.  Algorithm  7  and  Algorithm  7N  are  comparable,  but  the  former  takes  less 
times  to  obtain  a  solution.  The  run  time  of  the  CE  heuristics  is  smaller  for  case  D2 
than  for  C2,  see  Table  14.  This  is  counterintuitive  as  cases  D2  have  a  longer  time 
horizon  than  C2.  However,  the  randomness  in  the  algorithms  may  cause  such  effects. 
Overall,  it  appears  that  Algorithm  7  is  the  algorithm  of  choice  as  in  the  case  of  two 
searchers. 


Case 

Algo.  6 

CE  heuristic 
Algo.  7 

Algo.  7N 

Branch-and-bound  (Algo.  4) 
Optimal  Early  term. 

A2 

0.646586 

0.647400 

0.647400 

0.647400 

0.647400 

B2 

0.799565 

0.801566 

0.801552 

N/A 

N/A 

C2 

N/A 

0.336465 

0.336483 

N/A 

N/A 

D2 

N/A 

0.474782 

0.476186 

N/A 

N/A 

Table  13.  Probability  of  detection  for  problem  instances  with  two  searchers. 

(c)  Three  searehers 

Tables  15  and  16  report  the  numerical  results  for  the  cases  with  three 
searchers.  The  columns  are  the  same  as  the  corresponding  columns  in  Tables  11  and 
12.  We  observe  that  Algorithm  7  dominates  Algorithm  6  in  solution  quality  and  run 
times.  Furthermore,  it  is  available  for  the  large-scale  problem  instances.  As  in  the 
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Case 

Algo.  6 
(sec.) 

CE  heuristic 

Algo.  7  Algo.  7N 
(sec.)  (sec.) 

Branch-and-bound  (Algo.  4) 
Found  optimal  Proved  optimal 
(sec.) 

A2 

6.78 

11.09 

48.45 

0.67 

5.19 

B2 

35.23 

169.89 

517.05 

N/A 

N/A 

C2 

N/A 

650.93 

1666.43 

N/A 

N/A 

D2 

N/A 

210.43 

951.36 

N/A 

N/A 

Table  14.  Run  times  on  problem  instances  with  two  searchers. 


case  of  two  searchers,  Algorithms  7  and  7N  are  comparable,  but  the  former  tends  to 
require  less  computing  time  to  obtain  a  reasonable  solution. 

For  the  smallest  size  problem  instance  (A3),  Algorithm  7  obtains  an 
optimal  solution  (however  it  is  not  guaranteed)  in  49  seconds.  On  the  other  hand. 
Algorithm  4  obtains  an  optimal  solution  in  174.78  seconds  and  proves  it  optimal 
in  about  5  days.  The  best  detection  probability  after  49  seconds  in  Algorithm  4  is 
0.749631.  Thus  the  CE  heuristic  Algorithm  7  appears  to  generate  good  solutions 
quicker  than  the  B/B  based  Algorithm  4. 

Comparing  the  CE  heuristic  (Tables  11-16)  with  the  Static  Bound 
Heuristic  (SBH)  (Tables  8  and  9  second  to  last  row),  we  hnd  that  the  SBH  dom¬ 
inates  the  CE  heuristic  in  run  time  (3.25  seconds  compared  to  49  seconds  for  case 
A3).  However,  the  SBH  is  available  only  for  small-size  problem.  In  addition,  the  so¬ 
lution  quality  of  the  CE  heuristic  is  better  than  that  of  the  SBH  based  on  the  limited 
tests  . 

In  view  of  the  above  results.  Algorithm  7  appears  to  be  an  overall 
efficient  heuristic  for  solving  MSP. 
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Case 

Algo.  6 

CE  heuristic 
Algo.  7 

Algo.  7N 

Branch-and-bound  (Algo.  4) 
Optimal  Early  term. 

A3 

0.753688 

0.754400 

0.754400 

0.754400 

0.749631 

B3 

0.889145 

0.891720 

0.893104 

N/A 

N/A 

C3 

N/A 

0.436528 

0.436528 

N/A 

N/A 

D3 

N/A 

0.587159 

0.593178 

N/A 

N/A 

Table  15.  Probability  of  detection  for  problem  instances  with  three  searchers. 


Case 

CE  heuristic 

Algo.  6  Algo.  7  Algo.  7N 
(sec.)  (sec.)  (sec.) 

Branch-and-bound  (Algo.  4) 
Found  optimal  Proved  optimal 
(sec.) 

A3 

144.75 

49.00 

249.42 

174.78 

426,874.92 

B3 

767.39 

271.43 

1567.62 

N/A 

N/A 

C3 

N/A 

674.85 

3815.69 

N/A 

N/A 

D3 

N/A 

297.30 

1454.73 

N/A 

N/A 

Table  16.  Run  times  on  problem  instances  with  three  searchers. 
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IV.  MULTIPLE  HOMOGENEOUS 
SEARCHER  PROBLEM 

This  chapter  focuses  on  the  multiple-searcher  problems  where  all  searchers  are 
identical,  i.e.,  the  multiple  homogenous  searcher  problem  (MHSP).  The  goal  for  the 
searchers  is  to  minimize  the  probability  of  not  detecting  the  target  within  a  mission 
time,  which  is  equivalent  to  maximize  the  probability  of  detecting  the  target  within 
the  same  time.  We  utilize  the  convexity  of  the  nonlinear  objective  function  (the  non¬ 
detection  probability),  and  introduce  an  exact  algorithm  using  cutting  planes  (outer 
approximations).  Under  certain  assumptions,  the  problem  becomes  equivalent  to  a 
large-scale  linear  mixed-integer  program.  We  also  present  several  new  cuts  for  MHSP 
and  demonstrate  their  effect  in  computational  tests. 

A.  PROBLEM  DESCRIPTION  AND  FORMULATION 

Chapter  III  presented  three  algorithms  (the  branch-and-bound  procedure  and 
two  heuristics)  to  solve  the  multiple  searcher  problem  (MSP),  which  are  also  applica¬ 
ble  for  solving  MHSP.  Here  we  introduce  an  alternative  approach  especially  tailored 
to  solve  MHSP.  This  section  formulates  MHSP  by  following  [7,  18].  We  use  the  same 
notations  and  assumptions  as  in  MSP.  Specifically,  we  assume  that:  (1)  there  is  no 
end-point  restriction  for  any  searcher;  (2)  the  airspace  has  only  one  altitude,  i.e., 
there  is  only  one  vertex  corresponding  to  each  cell  in  the  discretized  area  of  interest; 
and  (3)  the  search  effect  at  the  current  waypoint  is  independent  of  the  previous  way- 
point.  In  addition,  we  assume  that  the  “search  effect”  is  described  by  an  exponential 
detection  function  instead  of  an  arbitrary  fnnction  as  in  earlier  chapters.  That  is, 
when  the  nnmber  of  searchers  in  cell  c  at  time  period  t  is  yc,t,  the  probability  of  de¬ 
tecting  the  target  in  that  cell  dnring  that  time  period  given  the  target  occupies  cell  c 
at  time  period  t  is  described  as  1  —  exp(— Q;c,ti/c,t)  instead  of  1  —  (1  —  g{c,  f))^'"’*  as  used 
in  Chapter  HI.  Here,  the  term  ac^t  is  the  detection  rate  for  a  single  searcher  in  cell  c 
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and  time  period  t  defined  by  ac^t  =  — log(l  —  g{c,  t)).  We  assume  that  g{c,  t)  <  1  (i.e., 
the  sensor  is  not  perfect)  so  that  ac,t  is  finite.  Additionally,  we  let  ui{t)  be  the  target’s 
cell  at  time  period  t  E  T,  the  sequence  of  cells  u  =  (ci;(l),  a;(2), . . . ,  ci;(T))  denote  a 
target  path,  and  be  the  probability  that  the  target  takes  that  path  u.  The  set  of 
all  possible  target  paths  is  denoted  as  hi.  In  this  notation,  MHSP  is  formulated  as 
follows. 

Formulation  of  MHSP 


Indices 

c,  c'  cells  (c,  c'  G  C  =  {1, . . . ,  C}) 
t  time  step  (f  G  T  =  {0, 1, ...,  T}) 

u  target  path  {u  G  hi) 

Parameters 

ac,t  detection  rate  for  a  single  searcher  in  cell  c  and  time  period  t 

ac,t  if  cell  c  is  on  target  path  u  at  time  period  t,  zero  otherwise. 
Pcfl  number  of  searchers  in  cell  c  at  time  period  0 
Pu)  probability  that  the  target  takes  path  u 

Variables 

Xc^c't  number  of  searchers  that  is  redistributed  from  cell  c  in  time 
period  t  to  cell  c'  in  time  period  t  +  I 
yc,t  number  of  searchers  in  cell  c  in  time  period  t 

Formulation 

min  f{y)  =  exp  -  ^ 

tJSO  \  c,t  J 

s.t. 

'y )  ^c',c,t—i  y 1) ...,T 

c'e7^(c)  c'e.F(c) 

yc,t  y  Xc^c',t  Vf  G  T 

c'e.5*='(c) 

Xc,c',t,  yep  integer 


We  refer  to  this  problem  which  minimizes  the  probability  that  the  target  is 
not  detected  during  time  horizon  T  subject  to  network  flow-balance  constraints  as 
the  multiple  homogeneous  searcher  problem  (MHSP).  The  objective  function  of  this 
nonlinear  mixed-integer  program  is  clearly  convex.  For  reference  later,  elements  of 
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the  gradient  V f{y)  of  the  objective  function  f{y)  are 

P<^ac,texp(- 

Ol/c,t  c,t 


(IV.l) 


Since  the  formulation  uses  all  possible  target  paths  ci;  G  ff,  if  the  number 
of  possible  target  paths  is  extremely  large,  calculating  the  objective  function  value 
/(y)  or  its  gradient  Vf{y)  for  any  solution  y  becomes  extremely  computationally 
expensive.  Brown  [7]  introduces  an  alternative  formulation  for  the  case  of  Markovian 
target  movement  using  conditioning  on  the  target’s  position  at  a  time  period  as 
follows. 


Given  a  solution  y,  let  Vc^tiy)  be  the  probability  of  the  target  visiting  cell  c  at 
time  period  t  without  being  detected  in  time  periods  1,  2,  ...,  t  —  1,  and  be  the 

probability  that  the  target  departs  cell  c  in  time  period  t  and  is  not  detected  by  any 
searches  in  time  periods  f  +  1,  t  +  2,  ...,  T.  Let  r.^tiy)  =  Vi,t{y)iT2,t{y)i  ■  ■  ■  )Tc,t{y)] 
and  s.^t{y)  =  S2,t{y),  •  •  • ,  sc,t{y)]-  We  dehne  r.^i{y)  =  p{-,  1),  where  p(-,  1)  is  a 

given  initial  target  distribution,  and  Sc^T{y)  =  1  for  any  cell  c  E  C.  Let  L  be  a  target 
transition  matrix.  Thus  r.^t{y)  and  s.^t{y)  are  calculated  recursively  by 


r.,t{y)  =  [ri,t-i{y)  exp{-ai^t-iyi,t-i),  •  •  • ,  rc,t-i{y)  exp{-ac,t-iyc,t-iW,  (IV.2) 

■S;tiy)  =  [spt+iiy)  exp(-ai,t+i|/i,t+i), . . . ,  rc,t+iiy)  exp{-ac,t+iyc,t+iW\  (IV.3) 

where  L'  is  the  transpose  matrix  of  L.  Since  we  only  consider  a  target  moving  ac¬ 
cording  to  a  Markov  transition  matrix,  we  can  apply  this  result.  In  this  notation,  for 
any  t  =  1,  2, ...,  T,  the  objective  function  is  alternatively  expressed  as 


fiy)  =  IlG,t(l/)exp(-ac,tl/c,t)sc,i(l/), 

cec 


and  elements  of  the  gradient  V  f{y)  are 

dfiy) 


dyc,t 


=  -OictrcAy)  exp(-ac,t|/c,t)sc,t(|/)- 


(IV.4) 


(IV.5) 


Thus,  for  any  solution  y,  the  objective  function  value  and  the  gradient  can  be  evalu¬ 
ated  with  a  moderate  computational  effort. 
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B.  ALGORITHMS  FOR  MHSP 

As  mentioned,  MHSP  is  a  convex  nonlinear  mixed-integer  program.  This  sec¬ 
tion  introduces  an  exact  outer  approximation  (OA)  algorithm  to  solve  MHSP  using 
the  convexity  of  the  objective  function  [30].  The  OA  algorithm  refers  to  the  fact  that 
the  surface  described  by  a  convex  function  lies  above  the  supporting  hyperplane  to 
the  convex  function  at  any  point.  A  supporting  hyperplane  at  a  point  takes  the 
form  /(yW)  +  V/(|/W)'(|/  -  yW)  (see  Figure  10).  The  algorithm  iteratively  generates 
supporting  hyperplanes  (cuts)  and  accumulates  them  to  provide  successively  improv¬ 
ing  linear  approximations  of  the  nonlinear  convex  function  (see  Figure  10).  The  linear 
approximation  underestimates  the  objective  function.  We  start  by  describing  a  basic 
OA  algorithm.  After  that  we  introduce  several  new  cuts  and  present  computational 
results  in  the  subsequent  subsections. 

1.  Basic  OA  Algorithm 

The  basic  OA  algorithm  is  described  as  Algorithm  8  (see  below),  which  solves 
the  following  master  problem  as  part  of  its  calculations.  Specihcally,  we  denote  that 
the  master  problem  of  the  /c-th  iteration  of  Algorithm  8  for  MP1(A;)  and  its  optimal 
value  and  optimal  solution  for  and  respectively. 

Formulation  of  Master  problem  :  MPl(/c) 


min 

s.t. 

Z  >  -  y^^)  i  =  0, 1, k-l 

c^G7^(c) 

yc,t  ^  ^  ^  ^ 

c'eJ^(c) 

Xc,c',t,  yc,t-  integer 

Algorithm  8  (Basic  OA  Algorithm). 

Step  0.  Set  a  relative  optimality  tolerance  5  >  0.  Choose  an  initial  feasible 
solution  y^^\ 
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Figure  10.  Linear  supporting  functions. 

Step  1.  Calculate  and  V/(|/‘^°^)  (see  (IV.4)  and  (IV. 5)).  Set  lower 

bound  q  =  0,  upper  bound  q  =  and  k  =  1. 

Step  2.  If  the  gap  {q—q)/q  <  S,  stop.  Else,  solve  the  master  problem  MP1(A;), 
and  obtain  its  optimal  value  and  optimal  solution  y^^\  If  z(^)  >  q,  then 
q  =  z^^\ 

Step  3.  Calculate  and  V f{y^^^)  (see  (IV.4)  and  (IV. 5)). 

Step  4.  If  <  q,  then  q  =  f{y^^^).  Replace  khj  k  +  1,  and  go  to  Step 

2. 

As  seen,  Algorithm  8  generates  one  cut  in  each  iteration  (a  strategy  we  refer 
to  as  “single-cut”)  at  a  solution  determined  by  the  master  problem  MP1(A;). 

It  is  well  known  that  using  a  nonzero  optimality  tolerance  when  solving,  at 
least  the  initial  master  problems  in  Algorithm  8  may  reduce  overall  computing  time 
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as  cuts  can  be  generated  more  cheaply.  We  examine  this  possibility.  Let  >  0  be 
a  relative  optimal  tolerance  for  the  solution  of  the  master  problem  MP1(A;).  Then, 

(k) 

the  value  obtained  from  the  master  problem  may  not  be  a  valid  lower  bound  on 
the  optimal  value  of  MHSP.  However,  (1  —  is  a  valid  lower  bound  and  we  use 

that  instead  of  in  step  2  of  Algorithm  8.  We  note  that  when  is  constant  at 
zero  for  iterations  after  some  hnite  iteration.  Algorithm  8  is  an  exact  algorithm  and 
provides  an  optimal  solution  after  a  hnite  number  of  iterations  [16]. 

We  examine  Algorithm  8  using  the  problem  instance  “case  C3”  in  Table  10 
of  Chapter  111  which  involves  15  by  15  cells,  a  time  horizon  of  18,  and  3  searchers. 
We  also  adopt  the  same  assumptions  and  parameter  settings  as  in  Chapter  111.  Our 
program  is  coded  using  the  General  Algebraic  Modeling  System  (GAMS)  Distribution 
22.5  [21]  and  is  implemented  on  the  same  computational  platform  as  in  Chapters  11 
and  111.  The  master  problem  MPl(/c)  is  solved  by  the  CPLEX  10.0  solver  [27].  In  the 
master  problem,  we  assume  that  the  optimality  tolerance  is  constant  with  value 
0  (version  1)  or  with  value  0.03  (version  2). 

All  searchers  depart  the  upper-left  corner  cell.  Thus,  as  a  choice  of  initial 
feasible  solution  (search  plan),  we  consider  the  situation  that  all  searchers  loiter  at 
the  depot  during  the  whole  time  horizon.  We  refer  to  this  search  plan  as  a  “trivial 
search”  plan.  This  trivial  plan  is  clearly  unwise,  but  we  apply  this  plan  as  an  initial 
feasible  solution  in  the  preliminary  tests.  Later  we  introduce  a  procedure  to  obtain  a 
better  initial  feasible  solution. 

In  the  initial  numerical  tests,  the  goal  is  to  compare  various  approaches  thus 
we  terminate  the  calculations  after  one  hour.  Table  17  reports  lower/upper  bounds 
on  the  optimal  value  and  the  relative  gap  between  the  bounds  at  every  ten  minute 
during  one  hour  calculations  using  Algorithm  8.  As  reference,  the  best  known  non¬ 
detection  probability  is  0.563472  (=1-0.436528)  which  is  found  in  Table  15  (case  C3) 
in  Chapter  111.  The  one  hour  run  time  includes  the  solution  time  for  the  master 
problems  (Step  2)  and  the  overhead  time  for  generating  cuts,  etc.  (Steps  1  and  3). 
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Algorithm  8  versions  1  and  2  use  96.4%  and  85.1%  of  the  run  time  to  solve  the  master 
problems,  respectively.  During  one  hour,  version  1  executes  only  25  iterations  because 
the  tighter  optimality  tolerance  in  the  master  problems,  while  version  2  executes  96 
iterations.  Even  though  the  cuts  in  version  2  may  not  be  as  “deep”  as  the  one 
generated  in  version  1,  it  is  clear  from  Table  17  that  the  extra  effort  to  solve  the 
master  problem  optimally  may  not  be  worth  it. 

In  these  test,  the  relative  optimal  tolerance  of  the  master  problem  is  con¬ 
stant  (0  or  0.03)  for  the  one  hour  calculations.  We  note  that  if  is  not  zero. 
Algorithm  8  may  regeneration  a  cut.  This  can  be  prevented  in  many  ways,  especially 
if  we  access  the  master  problem  solver.  However,  in  these  tests,  as  well  as  in  tests 
below,  we  adopt  the  simple  approach  of  adjusting  the  relative  optimality  tolerance 
Pqj^  example,  e^^^=0.03  appears  to  be  a  practical  number  if  the  goal  is  a  5% 
near-optimal  solution. 


Time 

Algo.  8 

version  1:  e 

(fc)  =  0 

Algo.  8  version  2: 

=  0.03 

(min) 

LB 

UB 

Gap 

LB 

UB 

Gap 

10 

0.410705 

0.620493 

0.511 

0.467784 

0.576832 

0.233 

20 

0.422249 

0.599865 

0.421 

0.471399 

0.576832 

0.224 

30 

0.428358 

0.599865 

0.400 

0.473591 

0.576832 

0.218 

40 

0.437659 

0.585815 

0.339 

0.475205 

0.574215 

0.208 

50 

0.448581 

0.585815 

0.306 

0.475316 

0.574215 

0.208 

60 

0.449400 

0.585815 

0.304 

0.475316 

0.572823 

0.205 

Table  17.  Numerical  results  for  Algorithm  8  versions  1  and  2  for  every  ten  minutes 
of  calculations.  The  best  known  non-detection  probability  is  0.563472  (=1-0.436528) 
found  in  Table  15  (case  C3)  in  Chapter  III. 


2.  New  Cuts  for  OA  Algorithm 

This  subsection  introduces  several  new  cuts  and  demonstrates  their  effect. 
Initially  two  new  cuts  (multi-cut  and  strong-cut)  are  presented.  After  that  we  prove 
that,  under  certain  assumption,  MHSP  is  equivalent  to  a  large-scale  linear  mixed- 
integer  program,  which  motivates  an  approach  for  obtaining  a  good  initial  solution. 
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Next,  two  additional  cuts  (relative-cut  and  symmetric-cut)  are  presented.  Finally,  we 
develop  a  specific  OA  algorithm  by  combining  these  cuts  effectively. 

(a)  Multi-cut 

In  Table  17,  we  learned  that  at  least  a  moderate  number  of  cuts  is  necessary 
to  bring  up  the  lower  bound.  The  aim  of  the  multi-cut  presented  here  is  that,  in 
each  iteration,  we  generate  multiple  cuts  instead  of  single  cut  to  accumulate  many 
cuts  fast  and  push  up  the  lower  bound  quickly.  The  multi-cut  utilizes  a  similar 
technique  to  the  multicut  version  of  the  L-shaped  method  in  stochastic  programming 
with  two-stage  recourse  (see  Chapter  5  in  [6]).  The  basic  idea  is  that,  for  the  objective 
function  f{y)  =  Puj fUv)  where  f^{y)  =  exp(- Ec,t  we  consider  outer 

approximation  of  foj{y)  for  all  a;  G  fl  instead  of  for  f{y),  and  generate  |f2|  cuts  at 
each  iteration.  However,  the  number  of  possible  target  paths  |f2|  is  extremely  large. 
Thus,  according  to  the  target  movement  during  the  initial  few  time  steps,  we  define 
U  (typically  U  «  |f2|)  “scenarios”  fli,  •••,  where  a  scenario  is  a  subset 
of  target  paths  (i.e.,  flu  (Z  fl,  u  =  1,2,...,  7/).  Moreover,  each  scenario  is  mutually 
exclusive  (i.e.,  fl  =  0,  m  7^  u')  and  each  possible  target  path  u  belongs  to 
some  scenario  (i.e.,  [J^^iflu  =  fl).  For  example,  from  a  known  target  location  at  time 
period  t  =  1  there  are  five  possible  target  movements  (stay  in  the  initial  position  or 
go  up /down/ left /right)  for  the  next  time  period.  Thus  we  can  define  five  {U  =  5) 
scenarios  based  on  the  target  movement  conditioned  during  the  initial  two  time  steps. 
Similarly,  for  the  initial  three  time  steps  (1=1,2  and  3),  twenty  five  scenario  {U  =  25) 
are  defined  if  boundary  effects  are  ignored.  Let  pu  be  the  probability  that  the  target 
takes  any  path  oj  such  that  u  &  flu- 

In  order  to  apply  this  multi-cut  approach,  the  objective  function  (IV. 4)  and 
the  gradient  (IV. 5)  need  to  be  expressed  differently.  When  the  solution  is  y,  let 
be  the  probability  of  target  visiting  cell  c  at  time  period  t  without  being  detected 
in  time  periods  1,  2,  ...,  1  —  1  given  the  target  takes  a  path  u  corresponding  to 
scenario  u.  Similarly,  s^Py)  is  defined  as  the  probability  that  the  target  departs  cell 
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c  in  time  period  t  and  is  not  detected  by  any  searches  in  time  periods  t  +  1,  t  +  2, 
T  given  the  target  chooses  a  path  uj  corresponding  to  scenario  u.  Let  r^t{y)  = 
,  r^,t{y)]  and  s'f^^{y)  =  [sltiy) ,  ^Itiy) ,  •  •  • ,  (?/)]•  %  considering 

the  conditioning  on  the  target’s  position  with  respect  to  scenario  u,  both  r^t{y)  and 
s^^{y)  are  recursively  calculated  similar  as  r.^t{y)  and  s.^t{y)-  In  this  notation,  the 
objective  function  is  redehned  as 

f{y)  =  '^Pufu{y),  (IV.6) 

u=l 


where,  for  any  t  =  1,2,  ...,T,  fu{y)  =  Ecec exp(-ac,t|/c,t)-s“t(|/),  and  elements 
of  the  gradient  V  fu{y)  are 


dfuiy) 

dyc,t 


exp(-ac,i|/c,tXt(l/)- 


(IV.7) 


The  OA  algorithm  with  U  cuts  per  iteration  is  described  as  Algorithm  9. 
The  algorithm  solves  the  following  master  problem  MP2(/c)  in  the  k-ih.  iteration. 
We  denote  the  corresponding  optimal  value  and  optimal  solution  for  and  y^^\ 
respectively. 

Formulation  of  Master  problem  :  MP2(/c) 


min  ^  =  Y!i=i  PuZu 
s.t. 

Zu  >  fuiy^"'^)  +  V/„(|/ (*))'(?/  -  !/(*))  u  =  l,2,...,U,  f  =  0, 1, ...,  /c  -  1 
5I/c'e7^(c)  —  X]c'eJP(c)  ^c,c',t  =  1, ...,  T 

yc,t  X/c'eJP(c)  ^c,c',t  ^  T 

Xc,c',t,  yc,t-  integer 


Algorithm  9  (OA  Algorithm  with  multi-cut). 

Step  0.  Set  a  relative  optimality  tolerance  5  >  0.  Choose  an  initial  feasible 
solution  y^^\ 

Step  1.  Calculate  fu{y^^^)  and  V/„(|/*'°^),  m  =  1, 2, ...,  U  (see  (IV.6)  and  (IV.7)). 
Set  lower  bound  q  =  0,  upper  bound  q  =  Y.u=iPufuiy^^^),  and  k  =  1. 
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Step  2.  If  the  gap  {q—q)/q  <  S,  stop.  Else,  solve  the  master  problem  MP2(A;), 
and  obtain  its  optimal  value  and  optimal  solution  If  zi^)  >  g,  then 
q  =  z^’^\ 

Step  3.  Calculate  and  V u  =  1,2, ...,  U  (see  (IV. 6)  and  (IV.7)). 

Step  4.  If  /(|/(^))  =  Y.u=iPufu{y^'''^)  <  g,  then  g  =  Replace  k  by 

k  +  1,  and  go  to  Step  2. 

Similar  to  the  master  problem  MP1(A;)  in  Algorithm  8,  the  relative  optimality 
tolerance  of  the  master  problem  MP2(/c)  is  denoted  by  >  0.  Again  we  note  that, 
if  >  0,  z^^^  must  be  replaced  by  (1  —  to  obtain  a  valid  lower  bound  in 

Step  2  of  Algorithm  9.  Here,  z^^})  is  the  value  turn  by  the  solver  when  using  relative 
optimality  tolerance  >  0  when  solving  MP2(A;). 

We  examine  Algorithm  9  using  the  same  problem  instance  as  in  Table  17.  In 
these  instances,  conditioning  on  the  target  movement  in  the  initial  two  time  steps 
{t=l  and  2)  results  in  U=5  scenarios.  Conditioning  on  the  initial  three  time  steps 
results  in  7/  =  25  scenarios.  We  refer  to  these  two  versions  of  Algorithm  9  as  “multi- 
t2-cut”  and  “multi-t3-cut,”  respectively.  As  above,  the  relative  optimality  tolerance 

when  solving  the  master  problem  MP2(/c)  is  set  to  0  or  0.03.  This  results  in  a 
total  of  four  versions  of  Algorithm  9: 

•  version  1:  multi-t2-cut,  and  MP2(A;)  solved  with  =  0 

•  version  2:  multi-t2-cut,  and  MP2(A;)  solved  with  =  0.03 

•  version  3:  multi-t3-cut,  and  MP2(A;)  solved  with  =  0 

•  version  4:  multi-t3-cut,  and  MP2(A;)  solved  with  =  0.03 

We  note  that  the  initial  feasible  solution  is  still  set  to  be  the  trivial  search 
plan  (i.e.,  all  searchers  keep  searching  the  initial  cell  for  the  whole  duration).  Table 
18  shows  the  numerical  results  for  versions  1  and  2  of  Algorithm  9.  Versions  1  and 
2  spend  96.3%  and  80.1%  of  the  run  time  to  solve  the  master  problems  (in  Step 
2),  respectively.  The  remaining  time  is  used  to  generate  cuts  and  build  models. 


86 


During  one  hour,  the  version  1  permits  only  20  iterations  while  version  2  executes 
68  iterations.  Since  the  multi-cut  requires  much  more  overhead  time  to  generate 
(multiple)  cuts  than  Algorithm  8  (single-cut)  in  each  iteration,  the  number  of  possible 
iterations  is  reduced  compared  to  that  of  Algorithm  8  (i.e.,  25  iterations  in  version  1 
of  Algorithm  8  are  reduced  to  20  in  version  1  of  Algorithm  9;  96  iterations  in  version 
2  of  Algorithm  8  are  reduced  to  68  in  version  2  of  Algorithm  9).  However,  for  both 
version  1  and  version  2  of  Algorithm  9,  the  total  nnmber  of  accnmulated  cuts  becomes 
larger  than  the  corresponding  versions  using  a  single-cnt  approach.  As  a  result,  the 
lower  bound  of  the  mnlti-t2-cut  (versions  1  and  2)  improves  quicker  than  that  of  the 
single-cnt  (refer  to  Table  17).  For  the  npper  bonnd,  there  is  no  signihcant  difference 
between  the  mnlti-t2-cut  and  the  single-cnt  approaches.  After  one  honr,  the  mnlti- 
t2-cnt  versions  provides  better  relative  gaps  than  the  single-cnt  versions  (i.e.,  30.4% 
(version  1  of  Algorithm  8)  compared  to  28.3%  (version  1  of  Algorithm  9),  and  20.5% 
(version  2  of  Algorithm  8)  compared  to  19.6%  (version  2  of  Algorithm  9)).  We  also 
note  that  version  2  =  0.03)  performs  better  than  the  version  1  (e^^^  =  0))  as  seen 

from  Table  17. 


Time 

Algo.  9  version  1:  e 

(fc)  =  0 

Algo.  9  version  2:  e^) 

=  0.03 

(min) 

LB 

UB 

Gap 

LB 

UB 

Gap 

10 

0.436285 

0.588137 

0.348 

0.470301 

0.569675 

0.211 

20 

0.454153 

0.588137 

0.295 

0.474848 

0.569675 

0.200 

30 

0.457303 

0.588137 

0.286 

0.475066 

0.569675 

0.199 

40 

0.458369 

0.588137 

0.283 

0.475835 

0.569675 

0.197 

50 

0.458369 

0.588137 

0.283 

0.476396 

0.569675 

0.196 

60 

0.458369 

0.588137 

0.283 

0.476396 

0.569675 

0.196 

Table  18.  Nnmerical  resnlts  for  Algorithm  9  using  mnlti-t2-cnt  (versions  1  and  2) 
for  every  10  minntes  of  calculations.  The  best  known  non-detection  probability  is 
0.563472  same  as  in  Table  17. 


Table  19  presents  the  nnmerical  results  for  Algorithm  9  using  mnlti-t3-cnt.  In 
each  iteration,  the  mnlti-t3-cut  reqnires  much  more  time  to  generate  25  cuts.  Thus 
the  versions  3  and  4  spend  only  40.1%  and  4.7%  of  the  rnn  time  on  solving  the 
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master  problems  (in  Step  2)  and  the  remaining  time  generating  cuts.  During  one 
hour  tests,  versions  3  and  4  permit  only  12  and  19  iterations,  respectively.  Typically, 
the  multi-t3-cut  versions  provide  a  signihcant  improvement  in  the  lower  bound  in  each 
iteration.  However,  the  long  run  time  per  iteration  prevents  them  from  carrying  out 
“enough”  iterations.  Thus  it  appears  that  the  multi-t3-cut  versions  are  inferior  to  the 
single-cut  and  multi-t2-cut  versions  (see  Tables  17  and  18).  We  note  however  that 
our  implementation  of  Algorithm  9  in  GAMS  could  be  improved  by  implementing  the 
time-consuming  cut  generation  in  C-|--|-  or  some  other  faster  programming  language. 
If  the  cut  generation  time  could  be  reduces,  the  multi-t3-cut  versions  may  prove 
competitive.  However,  in  this  dissertation,  we  do  not  consider  this  programming 
enhancement. 


Time 

Algo.  9  version  3:  e 

(fc)  =  0 

Algo.  9  version  4: 

=  0.03 

(min) 

LB 

UB 

Gap 

LB 

UB 

Gap 

10 

0.255538 

0.711572 

1.785 

0.247878 

0.716509 

1.891 

20 

0.359439 

0.624907 

0.739 

0.343276 

0.662075 

0.929 

30 

0.399356 

0.604164 

0.513 

0.384563 

0.608274 

0.582 

40 

0.424481 

0.604164 

0.423 

0.419195 

0.605206 

0.444 

50 

0.424481 

0.604164 

0.423 

0.435430 

0.589771 

0.354 

60 

0.429744 

0.594532 

0.383 

0.450207 

0.581592 

0.292 

Table  19.  Numerical  results  for  Algorithm  9  using  multi-t3-cut  (versions  3  and  4) 
for  every  10  minutes  of  calculations.  The  best  known  non-detection  probability  is 
0.563472  same  as  in  Table  17. 


Based  on  the  test  results  in  Tables  17,  18  and  19,  we  conclude  that  Algorithm 
8,  version  2  (single-cut  with  e^^^=0.03),  and  Algorithm  9,  version  4  (multi-t2-cut  with 
e*'^^=0.03),  are  the  most  promising.  Hence,  we  adopt  these  versions  as  baselines  for 
further  comparisons.  When  the  meaning  is  clear  from  the  context,  we  refer  to  those 
versions  simply  as  single-cut  and  multi-cut,  respectively. 

(b)  Strong-cut 

In  Algorithms  8  and  9,  at  each  iteration,  cuts  are  generated  at  the  optimal  so¬ 
lutions  y  obtained  from  the  master  problems.  The  cuts  are  the  tangent  hyper-planes 


y  9+1 


Figure  11.  Strong  cut  and  tangent  hyper-plane  on  an  exponential  objective  function 
in  mixed-integer  program. 

and  utilize  that  f(y)  >  f{y)  +  f{y)\y  —  y)  for  all  y.  (For  simplicity,  we  now  argue 
using  single-cut,  but  the  same  arguments  can  be  build  for  multi-cut  versions).  We 
now  use  the  fact  we  are  only  concern  with  integer  values  of  y  to  build  a  stronger  cut. 
Figure  11  illustrates  the  basic  idea  of  the  new  cut.  Since  the  MHSP  is  an  integer 
program,  we  utilize  hnite  differences  of  the  objective  function  f{y)  by  considering  the 
perturbation  from  t/c^t  to  yc,t  +  1  while  keep  all  other  variables  hxed.  Theorem  IV.  1 
presents  this  new  cut.  Let  Ac,t  G  be  a  vector  in  which  the  (c,  t)  element  is  one 
and  the  other  elements  are  all  zero,  where  Z  is  the  set  of  integers. 
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Theorem  IV.l  For  any  y,  y  e  ,  f{y)  >  /(y)  +  Ec,J/(y+ Ac,t) -/(y)](|/c,t-yc,t)- 


Proof.  f{y)  =  Ea.Pa.exp(-Ec,t<tl/c,t)  =  Ea;exp(-Ec,i  <tl/c,t+logPa.).  Let  be 
a  CT-dimensional  vector  defined  by  a‘^  =  (OcJ  and  =  (— logp^).  Then,  >  0  and 
b‘^  >  0.  Hence,  f{y)  =  fuiv)  where  =  exp(— a‘^1/ —  6*^).  In  this  notation,  the 
result  is  equivalent  to  f^{y)  +  Ec,J/a;(y  +  ^c,t)  -  fUy)]iyc,t  -  yc,t)  <  fM  for  ah  cu. 
Consequently,  /a;(y)[l  +  Ec,i(exp(-<,t)  “  -  ijcy)  -  exp{-a‘^{y  -  ^))]  <  0.  Let 

=  exp(— and  Af  be  the  set  of  the  cell-time  pair  (c,  t)  such  that  yc,t  —  yc,t  7^  0 
and  >  0  (i.e.,  cell  c  is  on  path  u  at  time  period  t).  Then  we  only  need  to  show 
that  EieAr(l  -  Pi){yi  -  m)  +  >  1,  where  i  =  (c,t)  G  A/”,  0  <  A  <  1. 

Let  ki  =  yi  -  jji,  then  we  can  define  ip{P)  =  EieAr(l  -  Pi)ki  +  Pi\  where 
P  =  (Pi)  is  a  |A/'|-dimensional  vector.  Let  p{P)  =  Then  Vp{P)  = 

(log(/5i), ...,  log(/5|;v'|))V(/^))  and  V‘^p{P)  =  A'Ap{P),  where  H  is  a  lA/”]  by  lA/”!  matrix 
in  which  the  first  low  is  (log(/5i), ...,  log(/3|Ar|))  and  the  other  elements  are  all  zero. 
Since  the  matrix  A' A  is  positive  semi-definite  and  p{P)  >  0,  p{P)  is  convex  on  the 
set  with  Pi  >  Q  for  all  i  G  Af.  Thus  v^(/3)  is  also  convex  on  the  set  with  /dj  >  0  for 
all  i  G  A/”.  When  Pi=l  for  Vi  G  A/”,  Vip{P)  =  {kp-l  +  PP^), ...,  k\j^\{-l  +  /5|]y^|))  =  0. 
Thus,  the  minimum  value  of  ^fiP)  is  1  when  Pi=l  for  Vi  G  Af.  Consequently,  we 
obtain  that  ^fiP)  >  1  which  proves  the  result.  □ 

We  refer  to  this  type  of  new  cut  as  ’strong-cut’.  The  calculation  of  finite 
differences  [f{y  +  Z.c,t)  —  f(y)]  is  as  easy  as  that  of  the  gradient  (IV.5).  Specifically, 

lf(y  +  f^c,t)  -  f(y)] 

=  rc',t(y)[exp(-ac',tyc',t  -  Oic,t)  -  exp{-ac' ,tyc' ,t)]sc' ,t{y) 

c'ec 

=  rc,t(|/)[exp(-Q;c,t(l/c,t  +  1))  -  exp{-ac,tyc,t)]sc,t{y)-  (IV.8) 

For  the  multi-cut  version,  for  each  scenario  u  =  the  finite  differences 

[fu{y  +  f^c,t)  —  fu{y)]  are  calculated  similarly: 
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(IV.9) 


[fu{y  +  Ac,t)  -  fu{y)] 

=  H  ^^c',t(l/)[exp(-ac',tl/c',i  -  ac,t)  -  exp{-ac',tyc',t)K',t{y) 

c'eC 

=  rlt{y)[exp{-ac,t{yc,t  +  1))  -  exp{-ac,tyc,t)]slt{y)- 

Algorithm  10  describes  a  single-cut  OA  algorithm  using  the  strong-cut.  Algo¬ 
rithm  10  is  a  modihcation  from  Algorithm  8  since  at  each  iteration  a  strong  cut  is 
generated  instead  of  a  tangent  hyper-plane.  In  the  k-th  iteration  of  Algorithm  10,  we 
solve  a  master  problem  MP3  (A;)  dehned  below,  where  the  optimal  value  and  optimal 
solution  are  denoted  and  y^’^\  respectively. 


Formulation  of  Master  problem  :  MP3(/c) 


min  z 

s.t. 

z  >  /(i/»)  +  UcAfiy^^  +  Ae,t)  -  /(i/»)](i/c,i  -  y^l)  *  =  0, 1, ...,  k  -  1 

5Zc'e7^(c)  ^c',c,t— 1  X/c'eJ^(c)  ^c,c',t  Vt  1,  ...,T 
yc,t  X/c'eJP(c)  ^c,c',t  Vt  G  T 
Xc,c',u  yc,t-  integer 


Algorithm  10  (OA  Algorithm  with  single,  strong-cut). 

Step  0.  Set  a  relative  optimality  tolerance  5  >  0.  Choose  an  initial  feasible 
solution  y^^\ 

Step  1.  Calculate  f{y^^^)  and  -|-  A^,*)  —  f{y^^^)  for  V(c,  t)  (see  (IV. 4) 

and  (IV. 8)).  Set  lower  bound  q  =  0,  upper  bound  q  =  and  A;  =  1. 

Step  2.  If  the  gap  {q—q)/q  <  S,  stop.  Else,  solve  the  master  problem  MP3(A;), 
and  obtain  its  optimal  value  z^^^  and  optimal  solution  y^^\  If  z^^'^  >  q,  then 
q  =  z^^\ 

Step  3.  Calculate  and  for  y{c,t)  (see  (IV. 4) 

and  (IV.8)). 

Step  4.  If  f{y^’^^)  <  g,  then  q  =  f{y^^^).  Replace  A;  by  A;  +  1,  and  go  to  Step 

2. 
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Algorithm  11  shows  the  multi-cut  version  of  Algorithm  10.  In  the  k-ih.  it¬ 
eration  of  Algorithm  11,  the  master  problem  MP4(A;)  dehned  below  is  solved.  The 
optimal  value  and  optimal  solution  of  MP4(A;)  are  denoted  and  respectively. 
In  each  iteration  of  Algorithm  11,  U  cuts  are  generated  at  once. 


Formulation  of  Master  problem  :  MP4(/c) 


min  z  =  Y!i=i  PuZu 
s.t. 

Zu  >  fu{y^^)  +  +  Ac,t)  -  fu{y^"'^)]{yc,t  -  yc}) 

n  =  1, 2, ...,  f/,  i  =  0, 1, ..., /c  —  1 

Ec'e7?.(c)  ^c',c,t— 1  Ec'eJP(c)  Vt  1,  ...,T 

yc,t  Ec'eJ^(c)  ^c,c',t  Vf  G  T 
Xc,c',t,  yc,t-  integer 


Algorithm  11  (OA  Algorithm  with  multi/strong-cut). 

Step  0.  Set  a  relative  optimality  tolerance  <5  >  0.  Choose  an  initial  feasible 
solution  y^^\ 

Step  1.  Calculate  and  +  Ac,t)  -  /«(|/^°^),  n  =  1,  2, V(c,  t) 

(see  (IV. 6)  and  (IV.9)).  Set  lower  bound  g  =  0,  upper  bound  q  =  En=i  Pufu{y^^^), 
and  k  =  1. 

Step  2.  If  the  gap  {q—q)/q  <  S,  stop.  Else,  solve  the  master  problem  MP4(A;), 
and  obtain  its  optimal  value  z^^^  and  optimal  solution  y^^\  If  z(^)  >  g,  then 
q  =  z^^\ 

Step  3.  Calculate /„(|/(^))  and /„(|/(^) -h  A^,*)  - /„(|/('')),  u  =  1,2,  ...,C,  V(c,t) 

(see  (IV.6)  and  (IV.9)). 

Step  4.  If  f{y^^'>)  =  Y!i=iPufu{y^'"'^)  <  q,  then  q  =  f{y^^^).  Replace  k  by 
k  +  1,  and  go  to  Step  2. 

We  demonstrate  the  effect  of  the  strong-cut  using  the  same  problem  instance 
accompanied  with  the  same  assumptions  and  parameter  settings  as  in  the  previous 
subsection.  The  trivial  search  plan  is  still  used  as  the  initial  feasible  solution  y^'^k 
The  relative  optimality  tolerance  of  the  master  problems  is  assumed  constant  at  0.03. 
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Table  20  describes  the  effect  of  the  strong-cut.  Compared  with  the  previous 
algorithms  (Algorithms  8  and  9),  Algorithms  10  and  11  implement  1.4  times  more 
iterations:  96  for  Algorithm  8  compared  to  135  for  Algorithm  10,  and  68  for  Algorithm 
9  compared  to  95  for  Algorithm  11.  Furthermore,  the  relative  gap  after  one  hour  is 
drastically  reduced  from  about  20%  to  10.5%  for  both  versions.  Since  the  increase 
of  the  number  of  iterations  and  the  progress  in  the  gap  are  essentially  identical  for 
both  versions,  we  conclude  that  the  strong-cut  is  quite  effective.  The  performance  of 
Algorithms  10  and  11  during  one  hour  of  calculations  is  almost  parallel  to  each  other. 


Time 

Algo.  10  :  single-cut 

Algo. 

11  :  multi-cut 

(min) 

LB 

UB 

Gap 

LB 

UB 

Gap 

10 

0.509594 

0.570631 

0.120 

0.510988 

0.571003 

0.117 

20 

0.511369 

0.570631 

0.116 

0.512539 

0.569664 

0.111 

30 

0.512452 

0.568363 

0.109 

0.513929 

0.569664 

0.108 

40 

0.512819 

0.568363 

0.108 

0.514792 

0.569664 

0.107 

50 

0.512847 

0.568363 

0.108 

0.515407 

0.569664 

0.105 

60 

0.514277 

0.568363 

0.105 

0.515407 

0.569664 

0.105 

Table  20.  Numerical  results  for  Algorithms  10  and  11  using  strong-cut  for  every  10 
minutes  of  calculations.  The  best  known  non-detection  probability  is  0.563472  same 
as  in  Table  17. 


(c)  Choice  of  Initial  Solution 

So  far,  as  a  choice  of  initial  feasible  solution,  we  have  used  the  trivial  search 
plan  (i.e.,  all  searchers  keep  loitering  at  the  initial  cell).  The  choice  of  the  initial  solu¬ 
tion  certainly  influences  the  OA  algorithms,  especially  during  the  initial  calculations. 
In  this  section,  we  aim  to  develop  a  “good”  initial  solution.  In  advance  of  this,  we 
hrst  introduce  a  theoretical  property  of  the  MHSP,  which  facilitates  the  calculation 
of  a  good  initial  solution. 

We  now  assume  that  the  detection  rate  is  time  and  space  independent,  i.e., 
a  =  ac,t  for  all  c  and  t.  Under  this  assumption,  we  show  that  the  MHSP  is  equivalent 
to  a  large-scale  linear  mixed-integer  program.  Let  =  J2c,t  ^ctyc,t  be  the  total 
effective  search  effort  given  search  plan  y  when  the  target  takes  path  u  &  Q.  Since 
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the  detection  rate  is  assumed  identical,  is  a  multiple  of  a  and  thus  =  ma  for 
some  m  =  0,1, JT.  This  is  a  similar  idea  to  [52],  The  maximum  possible  total 
effective  search  effort  is  the  number  of  searchers  ( J)  times  the  time  horizon  (T),  which 
occurs  when  all  searchers  take  the  path  ui.  In  this  notation,  we  claim  that  the  MHSP 
is  formulated  as  the  following  large-scale  linear  mixed-integer  program. 

Formulation  in  linear  mixed-integer  program:  LI 


Indices 

c,  c'  cells  (c,  c'  G  C  =  {1, . . . ,  C}) 
t  time  step  (f  G  T  =  {0, 1, ...,  T}) 

j  searchers  {j  E  J  =  {1, ...,  J}) 

oj  target  path  {u  G  ff) 

Parameters 

a  detection  rate  for  a  single  searcher  in  cell  c  and  time  period  t 

ac,t  if  cell  c  is  on  target  path  u  at  time  period  t,  zero  otherwise. 
i/cfl  number  of  searchers  in  cell  c  in  time  period  0 
probability  that  the  target  takes  path  u 

Variables 

Xc^c't  number  of  searchers  that  is  redistributed  from  cell  c  in  time 
period  t  to  cell  c'  in  time  period  t  +  1 
yc,t  number  of  searchers  in  cell  c  in  time  period  t 

Formulation 

min 

s.t. 

Zu)  >  +  m  —  me“"]  -|-  (1/Q;)e“'^"(e“"  —  l)Wui 

\/uj  G  ff,  m  =  0, 1...,  JT 
5Zc'e7^(c)  1  ^c'£j^{c)  ^C,c' ,t  'Jt  1,  ...,T 

yc,t  X/c'eJP(c)  'Jt  G  T 

Xc,c',t,  yc,t-  integer 

Theorem  IV.  2  If  the  detection  rate  is  identical  (i.e.,  a  =  ac,t  for  all  c  and  t),  the 
convex  nonlinear  mixed-integer  MHSP  and  the  linear  mixed-integer  program  LI  have 
identical  global  minimum  solutions. 

Proof.  Let  fu^iW^)  =  where  =  T.c,tC^f,tyc,t  and  f{y)  = 

only  takes  a  hnite  number  of  discrete  value  hL^  =  ma,  m  =  0,1, ...,  JT.  Thus 
fojiWoj)  can  be  described  as  a  piece-wise  linear  function  (see  Figure  12).  The  linear 
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Figure  12.  Piece-wise  linearlization  of  the  exponential  function  fcj{W^). 

approximation  between  ma  and  (m  -|-  1)q;  is  obtained  as  f(W^)  =  +  m  — 

me~°‘]  +  (1/Q;)e“™"(e““  —  1)W^  by  solving  the  linear  equation  {bo  +  6i[mQ;]  =  e"™" 
and  bo  +  bi[{m  +  l)a]  =  for  parameters  bo  and  bi.  □ 

The  linearization  of  MHSP  under  the  assumption  of  identical  detection  rates 
explicitly  depends  on  all  possible  target  paths.  Hence,  the  program  LI  may  become 
extremely  large.  Thus  we  cannot  directly  apply  this  formulation  to  solve  most  in¬ 
stances  of  MHSP.  However,  this  linearization  provides  useful  insight  to  obtain  a  good 
initial  solution. 

The  linearlization  of  MHSP  considers  the  total  effective  search  effort  IT^  on 
each  target  path  cu  G  H  individually.  The  alternative  approach  is  to  consider  the 
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total  effective  search  effort  aggregated  over  all  target  paths.  This  approach  allow  us 
to  utilize  the  Markovian  property  of  the  target  movement,  but  it  leads  to  a  relaxation 
as  we  now  described.  Recall  that  the  target  moves  between  cells  according  to  a 
known  transition  matrix  T.  Let  p{-,t)  =  t  >  1  be  the  undetected  target 

distribution  in  time  period  t,  (p(-,  1)  is  known).  Thus  the  aggregated  total  effective 
search  effort  is  described  as  W  =  J2c,tP{c,t)ayc,y  In  this  notation,  the  aggregation 
of  the  large-scale  linear  mixed-integer  program  LI  takes  the  form  of  the  moderately 
sized  linear  mixed-integer  program  L2: 

Formulation  in  linear  mixed-integer  program  (aggregation):  L2 


Indices 

c,  c'  cells  (c,  c'  G  C  =  (1, . . . ,  C}) 
t  time  step  (t  G  T  =  {0, 1, ...,  T}) 

j  searchers  {j  E  J  =  {1, ...,  J}) 

Parameters 

a  detection  rate  for  a  single  searcher  in  any  cell  and  time 
i/cfl  number  of  searchers  in  cell  c  in  time  period  0 
p{-,t)  undetected  target  distribution  at  time  period  t 

Variables 

Xc^c't  number  of  searchers  that  is  redistributed  from  cell  c  in  time 
period  t  to  cell  c'  in  time  period  t  +  1 
yc,t  number  of  searchers  in  cell  c  in  time  period  t 

Formulation 

min  z 

s.t. 

W  =  Ec,tP(c,  t)ayc,y 

z  >  +  m  —  me~°‘]  +  (l/a)e“'""(e“"  —  l)hL,  m  =  0, 1...,  JT 

Ec'&TZ{c)  ^c' X)c'eJ^(c)  Vt  1,  ...,T 

yc,t  Ec'GT{c)  Vt  G  T 

Xc,c',t,  Vcp  integer 

We  observe  that  solutions  of  L2  is  not  identical  to  solutions  of  LI  as  L2  is  a 
relaxation  of  LI  due  to  the  aggregation  of  all  the  target  paths.  However,  L2  is  easily 
solved  and  typically  provides  a  good  initial  solution  for  the  main  calculations.  We 
will  use  this  aggregated  program  L2  to  obtain  an  initial  solution. 
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In  addition,  we  also  consider  a  simple  condition  (constraint)  to  obtain  a  better 
initial  solntion.  If  the  number  of  searchers  is  moderate  (e.g.,  3  searchers)  and  the  size 
of  the  area  is  large,  it  is  typically  better  that  the  searchers  visit  cells  not  previously 
searched.  (Of  course,  the  effect  of  this  strategy  also  depends  on  the  movement  of  the 
target.)  Given  a  specihc  G  T  (e.g.,  t'=l  for  the  problem  instance  with  3  searchers), 
a  non-revisit  constraint  is  described  as  yc,t  <  1,  Vc  G  C.  Therefore,  in  the  initial 
step  (Step  0)  of  the  OA  algorithms,  we  solve  the  aggregation  problem  (L2)  with  this 
non-revisit  constraint  to  obtain  an  initial  solution.  This  initial  problem  is  referred  as 
L3  as  follows. 

Initial  problem:  L3 


Indices 

c,  c'  cells  (c,  c'  G  C  =  {1, . . . ,  G}) 
t  time  step  (t  G  T  =  {0, 1, ...,  T}) 

j  searchers  {j  E  J  =  {1, ...,  J}) 

Parameters 

a  detection  rate  for  a  single  searcher  in  any  cell  and  time 
Ucfi  number  of  searchers  in  cell  c  in  time  period  0 
p(-,t)  undetected  target  distribution  at  time  period  t 

Variables 

Xc,c't  number  of  searchers  that  is  redistributed  from  cell  c  in  time 
period  t  to  cell  d  in  time  period  t  +  I 
Uc^t  number  of  searchers  in  cell  c  in  time  period  t 

Formnlation 

min  z 

s.t. 

^  =  Ec,tP(c,  t)ayc,y 

z  >  +  m  —  me““]  -f  (1/Q;)e“'™"(e“"  —  1)1T,  m  =  0, 1...,  JT 

5I/c'e7^(c)  2^c',c,t-l  =  X]c'ej^(c)  dt  =  1, ...,  T 

yc,t  X)c'eJP(c)  dt  G  T 

Et>t'yc,t  <  1,  Vc  G  C 

Xc,c',t,  yc,t-  integer 


In  the  OA  algorithms,  the  hrst  cut(s)  is  generated  at  the  initial  solution  ob¬ 
tained  from  the  initial  problem  L3.  In  addition,  we  also  add  the  cut  at  the  solution 
y  =  0  (i.e.,  f{y)  >  /(O)  -I-  V/(0)'|/)  which  can  be  shown  to  relate  to  the  mean  bound 
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in  a  branch- and-bound  procedure  [34,  52],  The  cut  at  y  =  0  linearlizes  the  probability 
of  non-detection  at  the  point  of  “no  search”  effect  so  it  gives  an  estimate  of  initial 
effect  of  increasing  the  search  effort  above  zero  in  a  cell.  We  apply  the  strong-cut, 
not  the  gradient-based  hyper  plane,  at  y  =  0  and  the  initial  solution  found  when 
solving  L3.  We  note  that  this  initial  process  is  also  available  for  the  algorithm  using 
multi-cut  (Algorithm  11)  and  is  referred  as  the  initial  problem  L4  as  follows. 

Initial  problem:  L4 


Indices 

c,  d  cells  (c,  c'  G  C  =  {1, . . . ,  C}) 
t  time  step  (t  G  T  =  {0, 1, ...,  T}) 

j  searchers  {j  E  J  =  {1, ...,  J}) 

u  scenario  {u  =  1,  2, ...,  U) 

Parameters 

a  detection  rate  for  a  single  searcher  in  any  cell  and  time 
ycfl  number  of  searchers  in  cell  c  in  time  period  0 
pu  probability  that  the  target  takes  scenario  u 

undetected  target  distribution  at  time  period  t  when  scenario  u 

Variables 

Xc^c't  number  of  searchers  that  is  redistributed  from  cell  c  in  time 
period  t  to  cell  d  in  time  period  t  +  1 
Uc^t  number  of  searchers  in  cell  c  in  time  period  t 
Wu  total  effective  search  effort  when  scenario  u 

Formnlation 

min  YZ=iPuZu 

s.t. 

fW  =  Ec,tP“(c,t)a|/c,3/,  M  =  l,2,...,f/ 

Zn  >  e-™“[l  +  m  -  me-“]  +  (l/a)e-™“(e-“  -  1)W„ 

M  =  1,  2, ...,  f/,  m  =  0, 1...,  JT 
Z)c'e7^(c)  =  Hd&Tic)  Xc,c',t  dt  =  1, ...,  T 

yc,t  X)c'eJP(c)  dt  G  T 

Et>t’yc,t  <  1,  Vc  G  C 
Xc,d,t,  yc,t-  integer 

Table  21  reports  the  effect  of  the  choice  of  an  improved  initial  solution  using 
the  same  problem  instance  with  the  same  assumptions  and  parameter  settings  as 
above.  We  dehne  Algorithm  10.2  (single-cut)  and  Algorithm  11.2  (multi-cut)  when 
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an  initial  solution  is  computed  from  L3  in  Algorithm  10  and  from  L4  in  Algorithm 
11,  respectively.  For  both  Algorithms  10.2  and  11.2,  the  relative  gap  decreases  after 
10  minutes  from  0.107  to  0.103  for  Algorithm  10.2,  and  from  0.105  to  0.098  for 
Algorithm  11.2.  However,  after  30  minutes  of  calculations,  the  advantage  of  the 
improved  initial  solution  becomes  insignihcant  since  the  numbers  of  cuts  are  added 
and  the  cuts  largely  influence  the  solution  quality.  Algorithms  10.2  and  11.2  perform 
better  than  the  previous  cases  using  trivial  initial  solution,  thus  we  apply  them  as 
the  baseline  for  the  OA  algorithms  which  include  the  next  set  of  new  cuts. 


Time 

(min) 

Algo. 

LB 

10.2  :  single-cut 

UB  Gap 

Algo. 

LB 

11.2  :  multi-cut 

UB  Gap 

10 

0.510520 

0.568875 

0.114 

0.512064 

0.567808 

0.109 

20 

0.512370 

0.568875 

0.110 

0.513271 

0.567808 

0.106 

30 

0.512422 

0.568875 

0.110 

0.514057 

0.567808 

0.105 

40 

0.513546 

0.568875 

0.108 

0.514660 

0.567808 

0.103 

50 

0.513747 

0.568875 

0.107 

0.515650 

0.567808 

0.101 

60 

0.513747 

0.568875 

0.107 

0.515900 

0.567808 

0.101 

Table  21.  Numerical  results  for  Algorithms  10.2  and  11.2  using  an  improved  initial  so¬ 
lution  for  every  10  minutes  of  calculations.  The  best  known  non-detection  probability 
is  0.563472  same  as  in  Table  17. 


(d)  Relative- cut 

In  Table  21,  the  upper  bound  (i.e.,  the  best  feasible  value)  does  not  change 
after  ten  minutes.  Thus  it  is  possible  that  the  upper  bound  could  have  been  optimal. 
(Of  course,  it  is  not  in  this  case  as  we  know  a  better  value  from  Chapter  III).  Let  y 
be  the  best  solution  so  far  providing  the  upper  bound  q  =  f{y).  Since  the  objective 
function  is  convex,  y  is  guaranteed  to  be  optimal  if  f{y)  is  still  the  upper  bound  after 
cuts  are  added  at  all  neighbor  integer  points  y  oi  y  (the  Euclidian  distance  between 
y  and  y  is  one).  Thus,  we  explore  the  strategy  of  generating  cuts  around  the  best 
solution. 

Figure  13  illustrates  the  rule  that  determines  a  neighboring  integer  point  y 
where  we  generate  a  cut.  In  earlier  the  OA  algorithms,  each  iteration  obtains  a 
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Figure  13.  Illustration  of  cut  generation  around  the  best  solution. 

solution  y  and  generates  one  or  more  cuts  at  y.  Now,  we  also  generate  an  additional 
cut  at  y.  Thus,  in  the  algorithms  using  the  single-cut,  two  cuts  are  generated  in  each 
iteration.  The  upper  picture  in  Figure  13  shows  a  case  where  the  best  solution  so  far 
y  is  still  better  than  the  current  solution  y  (i.e.,  f{y)  <  f{y))-  Then,  we  determin 
an  integer  point  y  which  is  a  neighbor  of  y  and  is  the  closest  point  to  y  measured  in 
the  Euclidean  distance.  The  lower  picture  illustrates  the  opposite  case.  In  this  case, 
after  generating  the  cuts  at  y  and  y,  the  best  solution  y  is  updated  as  y.  We  denote 
the  cut  at  the  solution  y  as  a  “relative-cut.  ”  Moreover,  for  the  relative-cut,  we  also 
apply  the  strong-cut  instead  of  the  gradient-based  hyper-plane  cut. 

The  OA  algorithms  using  the  relative-cut  is  described  as  Algorithm  12  and 
Algorithm  13.  Algorithm  12  is  the  single-cut  version,  and  Algorithm  13  is  the  multi- 
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cut  version.  Algorithms  12  and  13  also  use  the  initial  problem  L3  and  L4  to  obtain  a 
reasonable  initial  solution,  respectively.  Thus  Algorithms  12  and  13  are  extension  of 
Algorithms  10.2  and  11.2  to  the  ones  with  relative-cut,  respectively. 

Algorithm  12  (OA  Algorithm  with  single/strong/relative-cut). 

Step  0.  Set  a  relative  optimality  tolerance  5  >  0  and  =  0.  Solve  the 
initial  problem  L3  {t'  =  1  for  the  problem  instance  with  3  searchers)  and 
set  the  solution  as 

Step  1.  Calculate  /(|/(°^),  f{y^^'>  +  Ac,t)  -  f{y^^'>)  and  f{y^^'^  +  /\c,t)  - 

for  V(c,t)  (see  (IV.4)  and  (IV. 8)).  Set  lower  bound  q  =  0,  upper 
bound  q  =  f{y^^^),  y  =  y^^\  and  k  =  2. 

Step  2.  If  the  gap  {q  —  q)/q  <  S,  stop.  Else,  solve  the  master  problem  MP3(/c) 
(described  above  Algorithm  10),  and  obtain  its  optimal  value  and 
optimal  solution  y  =  y^^\  If  >  g,  then  q  =  z^^\ 

Step  3.  Calculate  and  +  /2zc,t)  —  for  V(c,  t)  (see  (IV.4) 

and  (IV. 8)). 

Step  4.  Obtain  a  neighbor  point  y  around  the  best  solution  by  comparing 
f{y)  with  q  =  f{y).  Replace  k  hj  k  +  1.  Set  =  y. 

Step  5.  Calculate  and  +  /2zc,t)  —  for  V(c,  t). 

Step  6.  If  f{y)  <  q,  then  q  =  f{y)  and  y  =  y.  Replace  k  hy  k  +  1,  and  go  to 

Step  2. 

Algorithm  13  (OA  Algorithm  with  multi/strong/relative-cut). 

Step  0.  Set  a  relative  optimality  tolerance  5  >  0  and  =  0.  Solve  the 
initial  problem  L4  {t'  =  1  for  the  problem  instance  with  3  searchers)  and 
set  the  solution  as  y^^\ 

Step  1.  Calculate  fuiy^^'^),  fuiy^^'^),  Uiy^^^  +  ^c,t)  -  /«(l/^°^)  and  fuiy^^'^  + 

^c,t)  —  fuiy^^^),  for  u  =  1,2,...,!/,  V(c,  t)  (see  (IV.6)  and  (IV. 9)).  Set  lower 

bound  g  =  0,  upper  bound  g  =  Y.u=iPufu{y^^^),  y  =  y‘^^\  and  k  =  2. 

Step  2.  If  the  gap  (g  — g)/g  <  S,  stop.  Else,  solve  the  master  problem  MP4(/c) 
(described  above  Algorithm  11),  and  obtain  its  optimal  value  z^^'^  and 
optimal  solution  y  =  y^^\  If  z^^'^  >  g,  then  g  =  z^^k 
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step  3.  Calculate and +  A^,*) u  =  1,2,  ...,17,  V(c,t) 

(see  (IV.6)  and  (IV.9)).  Calculate  f{y)  =  T,^=iPufu{y) 

Step  4.  Obtain  a  neighbor  point  y  around  the  best  solution  by  comparing 
f{y)  with  q  =  f{y).  Replace  k  hj  k  +  1.  Set  y^'^^  =  y. 

Step  5.  Calculate and +  n  =  1,2,  ...,C,  V(c,t). 

Step  6.  If  fijj)  <  q,  then  q  =  f{y)  and  y  =  y.  Replace  k  hy  k  +  1,  and  go  to 
Step  2. 

We  test  the  effect  of  the  relative-cut  using  the  same  problem  instance  with 
the  same  assumptions  and  parameter  settings  as  above.  Table  22  describes  the  com¬ 
putational  results  of  Algorithms  12  and  13.  During  one  hour  of  calculations,  we  do 
not  identify  signihcant  effect  of  the  relative-cut  by  comparing  Table  21  and  Table  22. 
However,  the  relative-cut  strategy  seems  not  to  be  detrimental  and  it  allows  for  a 
faster  accumulation  of  cuts  which  may  be  benehcial  so  we  implement  the  relative-cut 
in  our  OA  algorithms. 


Time 

Algo. 

12:  single-cut 

Algo. 

13:  multi-cut 

(min) 

LB 

UB 

Gap 

LB 

UB 

Gap 

10 

0.510523 

0.571749 

0.120 

0.512066 

0.577300 

0.127 

20 

0.512079 

0.569864 

0.113 

0.514172 

0.570003 

0.109 

30 

0.513185 

0.569864 

0.110 

0.514721 

0.566784 

0.101 

40 

0.514212 

0.569864 

0.108 

0.514806 

0.566784 

0.101 

50 

0.514538 

0.569864 

0.108 

0.515639 

0.566784 

0.099 

60 

0.514538 

0.569864 

0.108 

0.515690 

0.566784 

0.099 

Table  22.  Numerical  results  for  Algorithms  12  and  13  using  relative-cut  for  every  10 
minutes  of  calculations.  The  best  known  non-detection  probability  is  0.563472  same 
as  in  Table  17. 


(e)  Symmetric-cut 

The  hnal  cut  we  present  is  quite  specihc  and  not  applicable  for  general  cases. 
However  it  is  somewhat  effective  in  specihc  situation.  Suppose  that  the  shape  of  the 
area  of  interest  (AOI)  is  symmetric  and  the  searchers  and  the  targets  are  initially  in 
the  cells  which  are  located  on  a  line  that  divides  the  AOI  symmetrically  (see  Figure 
14).  Suppose  also  that  the  detection  rate  is  identical  for  all  cells  and  time  periods. 
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And  finally,  suppose  that  the  target  moves  to  an  adjacent  cells  equally  likely.  We 
refer  to  this  specihc  situation  as  “symmetric.”  In  the  symmetric  situation,  for  a 
search  plan  V,  there  exist  a  symmetric  (or  mirror)  plan  V'  with  identical  probability 
of  nondetection,  i.e.,  q{V)  =  q{V'),  see  Figure  14. 

Thus,  the  basic  idea  of  the  symmetric-cut  is  that,  in  each  iteration,  cuts  are 
generated  not  only  at  the  current  solution  (plan)  V  but  also  at  its  symmetric  solution 
(plan)  v'  since  the  symmetric  solution  can  be  evaluated  from  the  current  solution. 
Algorithm  14  describe  the  OA  algorithm  using  the  symmetric-cut,  which  is  an  exten¬ 
sion  of  Algorithm  12.  In  Algorithm  12,  at  each  iteration,  two  cuts  are  generated  at 
the  current  solution  y  and  a  neighbor  point  y  around  the  best  solution  (by  comparing 
f{y)  with  /(!/)).  In  Algorithm  14,  two  cuts  are  additionally  generated  at  the  sym¬ 
metric  solution  y'  and  a  neighbor  point  y'  around  the  best  solution  (by  comparing 
fiy')  with  f{y)).  Thus,  in  Algorithm  14,  four  cuts  are  added  at  each  iteration. 

The  problem  instance  examined  in  numerical  test  so  far  is  symmetric.  Thus 
we  attempt  to  examine  the  effect  of  the  symmetric  cut  in  this  instance.  Since  the 
multi-cut  utilizes  the  conditioning  of  the  target  movement,  it  violates  the  symmetry 
property.  Thus  the  symmetric  cut  is  available  in  the  algorithm  using  the  single-cut 
only. 

Algorithm  14  (OA  Algorithm  with  single/strong/relative/symmetric-cut). 

Step  0.  Set  a  relative  optimality  tolerance  5  >  0  and  y^^'^  =  0.  Solve  the  initial 
problem  L3  (described  below  Algorithm  12,  and  t'  =  1  for  the  problem 
instance  with  3  searchers)  and  set  the  solution  as  y^^\ 

Step  1.  Calculate  /(|/(°^),  f{y^^'>  +  A^,*)  -  f{y^^'>)  and  f{y^^'^  +  /\c,t)  - 

for  V(c,t)  (see  (IV.4)  and  (IV. 8)).  Set  lower  bound  q  =  0,  upper 
bound  q  =  f{y^^^),  y  =  y^^\  and  k  =  2. 

Step  2.  If  the  gap  {q  —  q)/q  <  S,  stop.  Else,  solve  the  master  problem  MP3(/c) 
(described  above  Algorithm  10),  and  obtain  its  optimal  value  and 
optimal  solution  y  =  y^^\  If  >  g,  then  q  =  z^^\ 
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Search  plan  P 
Symmetric  plan  P 


Figure  14.  Symmetric  environment  where  the  symmetric-cut  is  available. 

Step  3.  Calculate  and  -|-  Ac,t)  —  for  V(c,  t)  (see  (IV. 4) 

and  (IV.8)). 

Step  4.  Obtain  a  neighbor  point  y  around  the  best  solution  by  comparing 
f{y)  with  q  =  f{y).  Replace  k  hj  k  +  1.  Set  =  y. 

Step  5.  Calculate  and  +  /^c,t)  —  f{y^^^)  for  V(c,  t)- 

Step  6.  For  the  solution  y,  obtain  the  symmetric  solution  y' .  Replace  k  by 

A;  -I-  1.  Set  y^^'^  =  y' . 

Step  7.  Calculate  and  +  /^c,t)  —  f{y^^^)  for  V(c,  t)- 

Step  8.  Obtain  a  neighbor  point  y'  around  the  best  solution  by  comparing 
f{y')  with  q  =  f{y).  Replace  khj  k  +  1.  Set  y^'^^  =  y' . 

Step  9.  Calculate  and  -|-  /\c,t)  —  f{y^^^)  for  V(c,  t)- 

Step  10.  If  f{y)  <  q,  then  q  =  f{y)  and  y  =  y.  Replace  khy  k  +  1  and  go  to 
Step  2. 
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Table  23  reports  the  effect  of  the  symmetric-cut  as  applied  to  the  same  problem 
instance  as  in  Table  22.  Comparing  with  Algorithm  12  in  Table  22,  the  relative  gap 
after  one  honr  of  calculations  is  reduced  from  0.097  to  0.092.  Both  the  lower  and 
upper  bound  are  improved.  Thus,  among  the  OA  algorithms  using  the  single-cnt. 
Algorithm  14  is  the  best  for  this  problem  instance. 


Time 

Algo.  14: 

Symmetric-cut 

(min) 

LB 

UB 

Gap 

10 

0.512580 

0.574371 

0.121 

20 

0.514023 

0.572793 

0.114 

30 

0.515390 

0.568628 

0.103 

40 

0.515390 

0.568628 

0.103 

50 

0.515812 

0.568628 

0.102 

60 

0.516447 

0.568628 

0.101 

Table  23.  Numerical  results  for  Algorithms  14  using  symmetric-cnt  for  every  10 
minntes  of  calculations.  The  best  known  non-detection  probability  is  0.563472  same 
as  in  Table  17. 


We  examined  several  OA  algorithms  in  which  the  new  cnts  (multi-,  strong-, 
relative-  and  symmetric-cnts)  are  cnmnlatively  applied.  Specihcally,  Algorithm  13 
(multi-cnt)  and  Algorithm  14  (single-cnt)  perform  best  among  the  algorithms  we 
examined.  For  the  problem  instance  (15  by  15  cells,  time  horizon  18  and  3  searchers) 
with  specific  assnmptions  and  parameters  settings  (e.g.,  identical  detection  rate  and 
symmetric  situation,  etc.).  Algorithms  13  and  14  provide  an  approximate  solntion 
with  relative  optimality  gap  about  10%  after  one  hour  of  calculations.  In  Chapter 
III,  for  the  same  problem  instance  (see  Tables  15  and  16),  the  Cross-Entropy  (CE) 
henristic  provides  a  better  solution  0.563472  (=1-0.436528)  after  674.85  seconds.  Thus 
the  CE  heuristic  performs  better  for  this  one-hour  computational  test,  bnt,  of  course 
the  CE  heuristic  does  not  provide  solution  quality  guarantees.  Thus  a  hybrid  method 
using  the  CE  henristic  and  the  OA  algorithm  is  expected  to  be  more  effective:  the 
CE  heuristic  is  initially  applied  to  obtain  a  good  initial  feasible  solntion  followed  by 
iterations  of  the  OA  algorithm.  This  dissertation  does  not  pursne  this  hybrid  scheme. 
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3.  Numerical  Study  of  Three  Searchers 

In  the  previous  subsection,  we  developed  specific  OA  algorithms  (i.e.,  Algo¬ 
rithm  13  and  Algorithm  14  among  others)  to  solve  the  MHSP.  For  the  moderately 
sized  problem  instance  (15  by  15  cells,  time  horizon  18  and  3  searchers)  examined, 
the  OA  algorithms  provided  a  solution  about  10%  of  the  optimal  solution  in  one  hour. 
We  now  examine  the  run  time  of  Algorithm  13  on  a  more  extensive  set  of  problem 
instances  of  somewhat  smaller  size.  Results  for  Algorithm  14  is  not  presented,  but 
they  are  almost  identical. 

As  problem  instances,  we  consider  the  following  five  cases  with  three  searchers 
in  which  we  apply  the  same  assumptions  and  parameter  settings  as  the  problem 
instance  in  the  previous  subsection,  except  for  changing  the  area  size  and  the  time 
horizon. 


•  case  1: 

•  case  2: 

•  case  3: 

•  case  4: 

•  case  5: 


15  by  15  cells  and  time  horizon  16 
13  by  13  cells  and  time  horizon  14 
11  by  11  cells  and  time  horizon  12 
9  by  9  cells  and  time  horizon  10 
7  by  7  cells  and  time  horizon  8 


For  each  case,  we  run  Algorithm  13  for  two  hours.  Since  the  problem  instances 
are  small,  the  master  problems  MP4(/c)  in  Algorithm  13  are  expected  to  be  solved 
quickly.  Thus,  we  allow  a  smaller  relative  optimality  tolerance  for  the  solver  of 
MP4(A;).  Specifically,  we  find  empirically  that  the  following  simple  rule  works  well; 
set  6*^^^=  min{0.03,(g  —  5')/(3g)}.  In  short,  at  each  iteration,  we  use  the  optimality 
tolerance  0.03,  but  apply  one  third  of  the  latest  gap  (g  —  q)/q  after  the  gap  has 
dropped  below  9%. 

Table  24  reports  the  gaps  at  every  ten  minute  during  two-hours  numerical  tests 
using  Algorithm  13  on  these  five  cases.  Essentially  all  cases  obtain  solution  within 
5%  of  optimality  within  20  minutes.  The  less  the  size  of  the  problem  instance  is. 
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the  smaller  gaps  are  achieved.  Case  5  results  in  an  essentially  optimal  solution  after 
two  hours.  After  two  hours,  the  lower  and  upper  bounds  in  that  case  are  0.454141 
and  0.455054,  respectively.  Thus  Algorithm  13  is  quite  effective  for  these  problem 
instances  and  provides  near-optimal  solutions  in  reasonable  computing  times. 


Gap 

easel  case2  case3  easel  ease5 


Time 

(min) 

To 

20 

30 

40 

50 

60 

70 

80 

90 

100 

no 

120 


0.056 

0.051 

0.049 

0.048 

0.048 

0.047 

0.046 

0.044 

0.041 

0.041 

0.041 

0.041 


0.050 

0.047 

0.044 

0.041 

0.040 

0.038 

0.038 

0.037 

0.037 

0.036 

0.036 

0.036 


0.039 

0.034 

0.031 

0.029 

0.028 

0.027 

0.026 

0.026 

0.026 

0.025 

0.025 

0.024 


0.025 

0.020 

0.018 

0.016 

0.015 

0.014 

0.013 

0.012 

0.012 

0.011 

0.010 

0.009 


0.022 

0.016 

0.013 

0.011 

0.009 

0.007 

0.006 

0.005 

0.004 

0.004 

0.003 

0.002 


Table  24.  Relative  optimality  gaps  for  Algorithm  13  applied  three  searcher  problem 
instances  during  two  hours  of  calculations. 


4.  Application  with  Large  Number  of  Searchers 

This  subsection  considers  the  MHSP  with  more  searchers,  especially  5,  10,  15, 
and  30  searchers.  Conceptually,  the  OA  algorithms  can  treat  the  HMSP  with  a  large 
number  of  searchers  more  easily  than  the  Cross-Entropy  (CE)  heuristics  since  they 
do  not  require  the  generation  of  specific  network  structures.  Moreover,  discussed  and 
demonstrated  below,  the  OA  algorithms  actually  may  beneht  from  more  searchers 
as  the  continuous  relaxation  of  the  master  problems  may  becomes  stronger  and  thus 
reduce  the  master  problem  solution  times.  In  this  subsection,  we  examine  the  effect 
on  the  OA  algorithms  with  respect  to  the  number  of  searchers.  For  the  computational 
tests,  we  use  the  same  problem  instance  (15  by  15  cells  and  time  horizon  18)  as  in 
the  previous  subsection  2,  except  we  change  the  number  of  searchers.  The  other 
assumptions  and  the  parameter  settings  are  identical  except  that:  (i)  The  identical 
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detection  rate  a  is  adjusted  according  to  the  number  of  searchers.  For  example,  for  the 
case  with  30  searchers,  detection  rate  O.la  is  used  to  make  the  search  effect  equivalent 
to  the  3  searcher  case  (i.e.,  3q;  =  30  ■  O.lo;).  (ii)  The  optimality  tolerance  for  the 
master  problem  solver  is  adjusted  adaptively  by  setting  min{0.03,(g  —  5')/(3g)}. 
For  the  cases  with  a  large  number  of  searchers,  the  master  problem  (mixed-integer 
program:  MIP)  becomes  almost  equivalent  to  the  continuous  relaxation  of  the  MIP. 
Thus,  the  master  problem  with  small  tolerance  (e.g.,  0.01)  can  be  solved  quickly  thus 
we  do  not  need  a  large  optimality  tolerance  (0.03).  On  the  other  hand,  the  MIP  with 
much  smaller  optimality  tolerance  (e.g.,  0.001)  is  quite  time-consuming  to  be  solved. 
(3)  In  the  initial  problem  to  hnd  a  good  initial  solution,  the  non-revisit  constraint  is 
ignored  since  that  constraint  may  be  detrimental  in  situations  with  a  large  numbers 
of  searchers  in  a  small  and  moderately  sized  areas. 

We  examine  the  following  four  cases  are  examined  on  the  problem  instance 
with  15  by  15  cells  and  time  horizon  18: 

•  case  1:  5  searchers 

•  case  2:  10  searchers 

•  case  3:  15  searchers 

•  case  4:  30  searchers 

Since  the  number  of  searchers  is  fairly  large,  a  continuous  relaxation  of  the 
MHSP  becomes  almost  equivalent  to  the  original  integer  MHSP  problem.  A  con¬ 
tinuous  relaxation  means  that  the  integrality  restriction  in  MHSP  is  relaxed,  which 
implies  that  the  searchers  can  be  “partitioned”  arbitrarily.  We  denote  the  instances 
corresponding  to  cases  1-4  with  the  continuous  relaxation  for  cases  lr-4r. 

For  the  computational  tests,  we  use  Algorithm  13  (multi-cut  version)  but 
Algorithm  14  generates  almost  identical  results.  We  refer  to  Algorithm  13  without 
integer  restrictions  in  the  master  problem  for  Algorithm  15.  Clearly,  Algorithm  15 
may  obtain  non-integer  solutions  and  thus  a  search  plan  that  is  not  implementable. 
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One  way  of  obtaining  an  (integer)  search  plan  from  the  non-integer  solution  is  by 
applying  some  rounding  heuristic.  However,  we  have  not  examined  that  approach  to 
obtain  integer  solutions  from  the  continuous  relaxation.  Algorithm  15  does  provide 
a  lower  bound  on  the  optimal  value  of  MHSP  as  it  solves  a  relaxation.  We  note  that 
gradient-descent  methods  may  also  solve  the  relaxed  MHSP.  However,  we  have  not 
pursued  that  avenue. 

Tables  25  to  28  report  the  results  of  one-hour  computational  test  for  the  case 
with  5,  10,  15  and  30  searchers,  respectively.  For  the  cases  with  more  than  5  searchers, 
after  one  hour.  Algorithm  13  provides  quite  good  solutions  with  gaps  at  about  2%. 
Even  for  5  searchers,  the  gap  is  moderate  about  5%.  Algorithm  15  (continuous  relax¬ 
ation)  gives  non-detection  probability  (NDP)  values  close  to  those  from  Algorithm 
13.  Especially  for  the  cases  with  large  number  of  searchers  (15,  30),  the  obtained 
NDP  values  from  Algorithms  13  and  15  are  almost  same. 

Comparing  with  the  lower  bounds  of  Algorithms  13  and  15,  we  hud  that 
Algorithm  15  provides  larger  lower  bounds  than  Algorithms  13.  Of  course,  the  upper 
bound  from  Algorithm  15  is  not  valid  for  MHSP.  Algorithm  15  solves  each  iteration 
more  quickly  than  Algorithms  13  since  it  solves  continuous  relaxation  problems.  Thus, 
Algorithm  15  generates  many  more  cuts  and  brings  up  the  lower  bound  quickly.  Thus, 
if  the  allowed  calculation  time  is  small,  a  hybrid  method  using  Algorithms  13  and 
15  may  be  better.  For  example.  Algorithm  15  (continuous  relaxation)  is  initially 
implemented  for  some  period  to  push  up  the  lower  bound  quickly  by  generating 
many  cuts,  and  later  Algorithm  13  is  implemented  to  provide  a  good  (integer)  search 
plan.  We  note  the  Algorithm  13  with  many  cuts  takes  much  run  time.  Thus,  this 
hybrid  approach  may  be  extremely  inefficient  when  Algorithm  15  has  generated  too 
many  “weak”  cuts  before  Algorithm  13  starts.  We  examined  this  hybrid  approach  for 
several  cases,  but  found  it  generally  to  be  an  inferior  approach  for  this  reason. 
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Time 

Algo.  13:  case 

1 

Algo. 

15:  case 

Ir 

(min) 

LB 

UB 

Gap 

LB 

UB 

Gap 

10 

0.505386 

0.554756 

0.098 

0.506761 

0.547174 

0.080 

20 

0.514396 

0.546065 

0.062 

0.518175 

0.538696 

0.040 

30 

0.516137 

0.546065 

0.058 

0.521317 

0.535791 

0.028 

40 

0.518325 

0.546065 

0.054 

0.523267 

0.534307 

0.021 

50 

0.518325 

0.546065 

0.054 

0.524436 

0.533078 

0.016 

60 

0.518810 

0.546065 

0.053 

0.525279 

0.532586 

0.014 

Table  25.  Numerical  results  for  Algorithms  13  and  15  for  the  case  with  5  searchers 
for  every  10  minutes  of  calculations. 


Time 

(min) 

Algo.  13:  case2 
LB  UB 

Gap 

Algo. 

LB 

15:  case2r 
UB 

Gap 

10 

0.501388 

0.549834 

0.097 

0.506429 

0.546805 

0.080 

20 

0.511533 

0.541721 

0.059 

0.518083 

0.538787 

0.040 

30 

0.518910 

0.535177 

0.031 

0.521340 

0.535598 

0.027 

40 

0.520400 

0.535177 

0.028 

0.523113 

0.534203 

0.021 

50 

0.521077 

0.535081 

0.027 

0.524351 

0.533240 

0.017 

60 

0.522048 

0.533785 

0.022 

0.525165 

0.532440 

0.014 

Table  26.  Numerical  results  for  Algorithms  13  and  15  for  the  case  with  10  searchers 
for  every  10  minutes  of  calculations. 


Time 

(min) 

Algo.  13:  case3 
LB  UB 

Gap 

Algo. 

LB 

15:  case3r 
UB 

Gap 

10 

0.495162 

0.553566 

0.118 

0.507577 

0.546699 

0.077 

20 

0.511244 

0.540650 

0.058 

0.517370 

0.538297 

0.040 

30 

0.517330 

0.537061 

0.038 

0.520760 

0.536154 

0.030 

40 

0.519869 

0.535204 

0.029 

0.522931 

0.534244 

0.022 

50 

0.521624 

0.534385 

0.024 

0.524083 

0.533416 

0.018 

60 

0.523123 

0.533276 

0.019 

0.524873 

0.532910 

0.015 

Table  27.  Numerical  results  for  Algorithms  13  and  15  for  the  case  with  15  searchers 
for  every  10  minutes  of  calculations. 
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Time  I  Algo.  13:  case4  I  Algo.  15:  case4r 


(min) 

LB 

UB 

Gap 

LB 

UB 

Gap 

10 

0.487626 

0.555553 

0.139 

0.508092 

0.542999 

0.069 

20 

0.505637 

0.542030 

0.072 

0.518134 

0.538112 

0.039 

30 

0.513595 

0.538593 

0.049 

0.521319 

0.535204 

0.027 

40 

0.516808 

0.537237 

0.040 

0.523258 

0.533536 

0.020 

50 

0.520377 

0.534938 

0.028 

0.524341 

0.533536 

0.018 

60 

0.522040 

0.533946 

0.023 

0.524980 

0.532239 

0.014 

Table  28.  Numerical  results  for  Algorithms  13  and  15  for  the  case  with  30  searchers 
for  every  10  minutes  of  calculations. 
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V.  CONCLUSIONS  AND  FUTURE 

RESEARCH 

A.  CONCLUSIONS 

This  dissertation  develops  models  and  solution  methodologies  to  solve  the 
discrete-time  path-optimization  problems  for  single  and  multiple  searchers  looking 
for  a  non-evading  moving  target  in  a  hnite  set  of  cells.  We  especially  focus  on  the 
following  problems: 

•  Single  searcher  problem  (SSP)  with  additional  resource  constraints  related  to 

risk  exposure  to  threats  and  fuel  consumption  (RSSP) 

•  SSP  for  multiple  searchers  (MSP) 

•  MSP  for  multiple  homogeneous  searchers  (MHSP) 

The  dissertation  starts  by  formulating  RSSP,  which  generalizes  existing  models  of  the 
SSP  by  considering  (i)  history-dependent  glimpse  detection  probability,  (ii)  multiple 
altitudes  for  the  searcher,  and  (iii)  multiple  constraints  on  “consumption”  of  resources 
such  as  time,  fuel,  and  risk. 

We  develop  a  specialized  branch-and-bound  (B/B)  algorithm  for  solving  RSSP 
and  propose  a  new  bound  (Lagrangian  directional  static  bound)  on  the  optimal  de¬ 
tection  probability  using  network  expansion  to  account  for  a  portion  of  the  history 
of  the  current  path  and  using  a  Lagrangian  relaxation  to  eliminate  resource  con¬ 
straints.  We  also  derive  a  series  of  network  reduction  procedures  that  tighten  the 
Lagrangian  relaxation  and  reduce  the  amount  of  path  enumeration  required  by  the 
B/B  algorithm. 

In  direct  comparison  with  a  state-of-the-art  algorithm  for  SSP,  the  proposed 
bound  and  network  reduction  procedures  reduce  the  run  times  by  at  least  one  order 
of  magnitude.  For  RSSP  with  time,  fuel,  and  risk  constraints,  as  well  as  two  altitudes, 
our  B/B  algorithm  solves  problem  instances  with  10  by  10  cells  and  a  time  horizon 
of  40  typically  within  20  minutes. 
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We  develop  a  B/B  algorithm  and  two  heuristics  (the  static  bound  heuristic 
and  the  cross-entropy  (CE)  heuristic)  to  solve  MSP.  Among  these  algorithms,  the  CE 
heuristic  performs  better  for  a  broad  range  of  problem  instances. 

We  also  focus  on  MHSP  and  develop  an  exact  outer  approximation  (OA) 
algorithms  using  several  novel  cutting  planes  (multi-cut,  strong-cut,  relative-cut,  and 
symmetric-cut).  In  empirical  studies,  this  algorithm  provides  search  plans  that  are 
guaranteed  to  be  within  5%  of  an  optimal  plan  in  less  than  about  20  minutes  for 
problem  instances  involving  15  by  15  cells,  3  searchers,  and  a  time  horizon  of  16. 
Instances  with  5  searchers  are  solved  even  faster.  In  addition,  we  prove  that  under 
certain  assumptions  the  nonlinear  multiple  homogenous  searcher  problem  is  equivalent 
to  a  large-scale  linear  mixed-integer  program. 

B.  FUTURE  RESEARCH 

The  CE  heuristic  appears  to  quickly  generate  a  good  feasible  solution  for  MSP. 
However,  it  does  not  guarantee  solution  quality.  A  possible  future  area  of  research 
would  be  to  develop  a  hybrid  algorithm  that  uses  the  CE  heuristic  to  generate  an 
initial  feasible  solution  for  an  OA  algorithm.  Currently  we  implement  the  CE  heuristic 
using  C-|--|-  and  execute  the  OA  algorithms  using  GAMS  with  the  CPLEX  solver.  The 
development  of  an  integrated  computational  program  is  worthy  to  effectively  solve  the 
MSP/MHSP.  Eurthermore,  an  improved  implementation  of  the  OA  algorithms  outside 
GAMS  may  reduce  cut  generation  time  and  prove  worthwhile. 

This  dissertation  does  not  consider  the  multiple  searcher  problem  with  resource 
constraints  (MRSP).  In  practical  applications,  the  mission  planner  would  like  to  hnd 
an  optimal  search  plan  for  multiple  “resource-constrained”  searchers.  A  combined 
approach  of  the  RSSP  algorithm  (Lagrangian  directional  static  bound)  with  the  MSP 
algorithm  (CE  heuristic)  could  be  developed  for  this  purpose.  Eor  example,  such 
an  approach  is  applicable  for  the  case  where  the  multiple  searchers  start  the  search 
operation  all  at  the  same  time  and  hnish  their  mission  at  the  same  time. 
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